{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "bd9929be",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "%pip install \"pybamm[plot,cite]\" -q    # install PyBaMM if it is not installed\n",
    "# import dependencies\n",
    "import matplotlib.pylab as plt\n",
    "import numpy as np\n",
    "import scipy.optimize\n",
    "\n",
    "import pybamm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1223d98",
   "metadata": {},
   "source": [
    "# Sensitivities and data fitting using PyBaMM\n",
    "\n",
    "PyBaMM input parameters [`pybamm.InputParameter`](https://docs.pybamm.org/en/stable/source/api/expression_tree/input_parameter.html) can be used to run many simulations with varying parameters. Here we will demonstrate PyBaMM's ability to calculate the senstivities of model outputs with respect to input parameters. \n",
    "\n",
    "To be more specific, given a model output $f(a)$, where $a$ is an input parameter, we wish to calculate $\\frac{\\partial f}{\\partial a}(a)$.\n",
    "\n",
    "First we will demonstrate using a toy model, given by the equations\n",
    "\n",
    "$$\\frac{dy}{dt} = a y$$\n",
    "\n",
    "with a scalar state variable $y$ and a scalar parameter $a$, and initial conditions\n",
    "\n",
    "$$y(0) = 1$$\n",
    "\n",
    "We will also define a model output given by $f = y^2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "25970cdf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# setup a simple test model\n",
    "model = pybamm.BaseModel(\"name\")\n",
    "y = pybamm.Variable(\"y\")\n",
    "a = pybamm.InputParameter(\"a\")\n",
    "model.rhs = {y: a * y}\n",
    "model.initial_conditions = {y: 1}\n",
    "model.variables = {\"y squared\": y**2}\n",
    "\n",
    "solver = pybamm.IDAKLUSolver(rtol=1e-10, atol=1e-10)\n",
    "t_eval = np.linspace(0, 1, 80)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f3d61bf",
   "metadata": {},
   "source": [
    "Note that we have used the [`pybamm.IDAKLUSolver`](https://docs.pybamm.org/en/stable/source/api/solvers/idaklu_solver.html) solver for this example, this is currently the recommended solver for calculating sensitivities in PyBaMM.\n",
    "\n",
    "We can solve the model using a specific value of $a = 1$. However, now we will also calculate the forward sensitivities of the model by setting the argument `calculate_sensitivies=True`. Note that this argument will also accept a list of input parameters to calculate the sensitivities for, but setting it to `True` will calculate the sensitivities for **all** input parameters of the model "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "1b1781a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "solution = solver.solve(model, [0, 1], inputs={\"a\": 1}, calculate_sensitivities=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d6f176e",
   "metadata": {},
   "source": [
    "We can now access the solution as normal, and the sensitivities using the syntax: `solution[output_name].sensitivities[input_parameter_name]`\n",
    "\n",
    "Note that a current restriction to the sensitivity calculation is that it will only return the sensitivities at the values of `t_eval` used to solve the model. Any interpolation between these values will have to be done manually"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "bf0a2d9c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnUAAAHWCAYAAAARl3+JAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABZNUlEQVR4nO3deViU5f4/8PcMMMM+CLIKCIIriuKaS6lpLqnZnmVmtlqe0+KpY55TlqdFW34dO+XR6phaarZ81WxRU3MtVxQFFwRBZBEQlBnWYZi5f38MTJGorPM888z7dV1zFcMz8n4cuflwryohhAAREREROTS11AGIiIiIqOVY1BEREREpAIs6IiIiIgVgUUdERESkACzqiIiIiBSARR0RERGRArCoIyIiIlIAFnVERERECuAqdYC2ZLFYkJeXBx8fH6hUKqnjEFEjCSFQWlqKsLAwqNX83bOp2PYROa6WtH+KLury8vIQEREhdQwiaqbs7GyEh4dLHcPhsO0jcnzNaf8UXdT5+PgAsP7F+Pr6SpyGiBrLYDAgIiLC9j1MTcO2j8hxtaT9U3RRVzfs4Ovry4aNyAFx6LB52PYROb7mtH+crEJERESkACzqiIiIiBSARR0RERGRArCoIyIiIlIAFnVERERECsCijoiIiEgBWNQRERERKQCLOiIiIiIFYFFHREREpAAs6oiIiIgUgEUdERERkQKwqCMiIiJSABZ1RERERArAoo6IiIhIAVjUEZFdfH0oG2/+eBIpuXqpoxAR2dXcdcexeEc69BWmNv06rm36pxMR1Vp3NAf7My4hqr0XenbQSR2HiMguisuM+PJgNgDgocEd2/RrsaeOiNpclcmMI1klAIAhMe2lDUNEZEfJtaMTnQK94OPu1qZfi0UdEbW5xKzLqDZbEKpzR1SAp9RxiIjsJjnHWtT1ssMIBYs6Impzv50tAgAMjgmASqWSOA0Rkf3U9dSxqCMiRfjtbDEADr0SkfOpK+riw/3a/GuxqCOiNlVaZcLx2uGHwTEBEqchIrKfi6VGXNBXQaUC4sJ82/zrsagjojZ16NwlmC0CUQGe6ODnIXUcIiK7qdvCKSbQG17att9whEUdEbWp39KtQ6/spSMiZ1M3ShFvp22cWNQRUZuqm083mPPpiMjJJOeWAIDd9uZkUUdEbeZyeTVOXjAAAAZ3Yk8dETmX3xdJsKgjIge3P8PaS9cl2BuBPlqJ0xAR2U+hoQoFBiPUKqCHHRZJACzqiKgN/Vq7Px23MiEiZ1PXSxcb5A1PjX1OZWVRR0Rt5tfaRRJDY1nUEZFzOW47ScLPbl+TRR0RtYnckkpkFpXDRa3CoE7+UschIrKr30+SsM/QK8CijojayK/p1qHX+HAdfNv4EGsiIjkRQvxe1NnhJIk6LOqIqE3UFXXDOPR6hd27d2PSpEkICwuDSqXChg0brnrtzJkzoVKpsGjRIrvlI6KWKTAYcbHUCBe1Cj1C2VNHRA5MCMH5dNdQXl6O3r17Y/Hixde8bv369di/fz/CwsLslIyIWsPxnBIAQOcgb3hoXOz2de2zHIOInMqZgjIUlRnh7qZGQqSf1HFkZ/z48Rg/fvw1r8nNzcVf//pXbNmyBRMmTLBTMiJqDSm2+XT22Z+uDnvqiKjV7a0deh0YHQCtq/1+S1UKi8WCadOm4cUXX0RcXJzUcYioiY7bedPhOrIt6qKioqBSqa54zJo1S+poRHQdv8+n4ykSzfH222/D1dUVzzzzTKOuNxqNMBgM9R5EJA0hhK2nzl7Hg9WR7fDroUOHYDabbR+npKTglltuwT333CNhKiK6HpPZggMZnE/XXImJifjggw9w5MgRqFSqRr1mwYIFmD9/fhsnI6LGuKCvQlFZNVzVKnS34yIJQMY9dYGBgQgJCbE9fvjhB8TExGD48OFSRyOia0jKLkF5tRntPN3QPcS+DZoS7NmzB4WFhYiMjISrqytcXV2RlZWFv/3tb4iKimrwNXPnzoVer7c9srOz7RuaiGzqNh3uEuwDdzf7Tj+RbU/dH1VXV2PVqlWYPXv2NX9zNRqNMBqNto85BEFkf3vSrEOvQ2PbQ61uXE8T/W7atGkYPXp0vefGjh2LadOmYcaMGQ2+RqvVQqvl2bpEcpCcWwLA/oskAAcp6jZs2ICSkhI8/PDD17yOQxBE0tubdhEAcGNnDr1eTVlZGdLT020fZ2ZmIikpCf7+/oiMjERAQP25iG5ubggJCUHXrl3tHZWImigpuwQA0DvCz+5fW7bDr3+0bNkyjB8//rp7NXEIgkha+kqTrUEb1jlQ2jAydvjwYSQkJCAhIQEAMHv2bCQkJGDevHkSJyOilrBYBI5nW4dfe0ewp+4KWVlZ2LZtG9atW3fdazkEQSStfWeLYRFAp0AvdPDzkDqObI0YMQJCiEZff+7cubYLQ0StJqOoDKXGGni4uaBrsI/dv77se+qWL1+OoKAgbr5J5AD21A693sReOiJyQknZv2867Opi/xJL1kWdxWLB8uXLMX36dLi6yr5Tkcjp7eV5r0TkxJKyLwOQZugVkHlRt23bNpw/fx6PPPKI1FGI6DrOF1cgq7gCrmoVbojhpsNE5HyO1fbU9YloJ8nXl3X315gxY5o074SIpLMn3Tr02jeyHby1sm5aiIhaXZXJjFMXrFupsaeOiBza3tr96YZxKxMickIn8vSosQi099ZKtlCMRR0RtViN2WI775X70xGRM0qyDb3qGn3EX2tjUUdELXYsRw9DVQ10Hm6ID/eTOg4Rkd3V7dHZR4JNh+uwqCOiFtt9xjqfblhse7jwaDAickLHJDxJog6LOiJqsV21Rd1NXTj0SkTOp7jMiPOXKgBA0tEKFnVE1CIlFdU4nlMCALipCzcdJiLnczzHOp8uJtALOg83yXKwqCOiFtmbXgSLALoEeyNUx6PBiMj5HJXB0CvAoo6IWqhuPh2PBiMiZ1W3SCKBRR0ROSohBHafsW5lwqFXInJGQghZLJIAWNQRUQucKShDvqEKWlc1Bkb7Sx2HiMjuzhVXQF9pgsZVjW4hvpJmYVFHRM1WN/R6Q6cAuLu5SJyGiMj+6nrp4sJ8oXGVtqxiUUdEzbY7rW4rEw69EpFzksOmw3VY1BFRs1RU1+BAxiUAwHAWdUTkpFjUEZHD23e2GNVmC8LbeSAm0EvqOEREdmesMeNkngEAizoicmA7U61DryO6Bkp2eDURkZRO5BlQbbagnacbIv09pY7Doo6Imk4IgZ1nCgEAI7oESZyGiEgaR7IuAwD6dWwni19uWdQRUZNlFJUj+1IlNC5qDIkNkDoOEZEkjpy3FnV9O7aTOIkVizoiarK6odeB0f7w1LhKnIaIyP6EEEis66mLZFFHRA5qZ2rt0GtXrnolIueUW1KJAoMRrmoV4sP9pI4DgEUdETVRZbUZBzKtW5mwqCMiZ3XkfAkAoEeYLzw08th8nUUdETXJvowiVNdY0MHPAzGB3lLHISKSRN0iib4yGXoFWNQRURPtOM2tTIiI6ubTyWWRBMCijoiaQAiBX05b59Pd3I1bmRCRc6qorsHJC9ZNh/uxqCMiR5RWWIbckkpoXdUYEtNe6jhERJI4nqOH2SIQ4uuOMJ271HFsWNQRUaPV9dINjgmQzcRgIiJ7S5TZpsN1WNQRUaNx6JWI6PdFEgmRftIG+RMWdUTUKPoKk+2305FdWdQRkXMSQthOkpDTfDqARR0RNdLutIswWwRig7wRIYODq4mIpJBZVI7LFSZoXNWIC9NJHaceFnVE1Cg7OPRKRGQbsYjvoIPGVV5llLzSEJEsmS0CO89Y96fj0CsROTO5Dr0CLOqIqBGO5ZTgUnk1fNxd0T9Kfg0ZEZG9HMkqASCvTYfrsKgjouv65ZR16PWmLoFwc2GzQUTOSV9pwpnCUgDyOh6sDltnIrqubacKAACju3PolYicV1J2CYQAIv09EeijlTrOFVjUEdE15VyuwOn8UqhVwIguLOqIyHn9cdNhOWJRR0TXVLfhcL+O7dDOSyNxGiIi6SRmXQIA9JXZpsN1WNQR0TVtq51PN6p7sMRJlGP37t2YNGkSwsLCoFKpsGHDBtvnTCYT5syZg169esHLywthYWF46KGHkJeXJ11gIoLJbLEtkhgQ7S9tmKtgUUdEV1VmrMH+s8UAOJ+uNZWXl6N3795YvHjxFZ+rqKjAkSNH8Morr+DIkSNYt24dUlNTcdttt0mQlIjqnMgzoNJkhs7DDV2CfKSO0yBXqQMQkXztTStCtdmCjgGeiAn0ljqOYowfPx7jx49v8HM6nQ5bt26t99xHH32EgQMH4vz584iMjLRHRCL6k0OZ1qHXAVHtoFarJE7TMPbUEdFVba9d9TqqWzBUKnk2Ys5Ar9dDpVLBz89P6ihETuuAraiT59ArwJ46IroKs0VgR6p1Ph2HXqVTVVWFOXPm4P7774evr2+D1xiNRhiNRtvHBoPBXvGInILFInC4dpGEXOfTAeypI6KrSMq+jKKyulMk5NuIKZnJZMK9994LIQSWLFly1esWLFgAnU5ne0RERNgxJZHypV8sQ0mFCe5uavQM00kd56pY1BFRg7aetPbSjewaJLtDq51BXUGXlZWFrVu3XrWXDgDmzp0LvV5ve2RnZ9sxKZHyHawdek2IaCfr9pDDr0TUoK0n8wEAt/TgVib2VlfQpaWlYceOHQgICLjm9VqtFlqt/Ha3J1KKQ+fkP/QKsKgjogZkXCzD2YvlcHNRYXjXQKnjKE5ZWRnS09NtH2dmZiIpKQn+/v4IDQ3F3XffjSNHjuCHH36A2WxGfr61wPb394dGww2gieytbuXrIBZ1RORotp60rnq9oVMAfN3dJE6jPIcPH8bIkSNtH8+ePRsAMH36dLz22mvYuHEjAKBPnz71Xrdjxw6MGDHCXjGJCNajEvP0VXBVq5Ag05Mk6rCoI6IrbKvdyoRDr21jxIgREEJc9fPX+hwR2Vfd0GtcBx08NfIum+Q72w9Abm4uHnzwQQQEBMDDwwO9evXC4cOHpY5FpGjFZUbbodWjeTQYETm5g5nW9nBgVDuJk1yfbEvOy5cvY+jQoRg5ciQ2bdqEwMBApKWloV07+f+lEjmy7acLYRFAzw6+CPPzkDoOEZGkDmZaj0qU86bDdWRb1L399tuIiIjA8uXLbc9FR0dLmIjIOdTNp2MvHRE5u+IyI85eLAfgGEWdbIdfN27ciP79++Oee+5BUFAQEhIS8Omnn0odi0jRKqvN2JN2EQAwpkeIxGmIiKR16Jx16LVzkDfaecl/5blsi7qMjAwsWbIEnTt3xpYtW/DUU0/hmWeewcqVK6/6GqPRCIPBUO9BRI2368xFVJksiPD3QPdQH6njEBFJylH2p6sj2+FXi8WC/v3746233gIAJCQkICUlBUuXLsX06dMbfM2CBQswf/58e8YkUpSfT1j3QxvTIwQqlUriNERE0qor6uS+P10d2fbUhYaGokePHvWe6969O86fP3/V1/CoHKLmM5kttq1MxsZx6JWInFu5sQYn8qwjfo4wnw6QcU/d0KFDkZqaWu+5M2fOoGPHjld9DY/KIWq+AxmXYKiqQYCXBv06cpU5ETm3xKzLMFsEOvh5OMxOALLtqXv++eexf/9+vPXWW0hPT8eaNWvwySefYNasWVJHI1KkLSd+P+vVRc2hVyJybvsyrFuZ3NDp2mcvy4lsi7oBAwZg/fr1+PLLL9GzZ0+8/vrrWLRoEaZOnSp1NCLFsViEbSuTMXHcyoSIaL+tqHOMoVdAxsOvADBx4kRMnDhR6hhEinc8V498QxW8NC4YEtNe6jhERJIqM9bgeI4eADA4hj11RORA6oZeR3QLgrubi8RpiIikdejcJZgtAhH+Hghv5yl1nEZjUUfk5IQQ2JxiLerGcdUrEdHvQ6/RjtNLB7CoI3J6qQWlyCwqh8ZVjZHdgqSOQ0Qkuf1nrUWdIw29AizqiJzepmRrL91NnQPhrZX1NFsiojZnqDIhOdc6n86RVr4CLOqInF7d0Ov4nhx6JSI6lHkJFgF0DPB0mP3p6rCoI3JiGRfLkFpQCle1CqO7cysTIqK6+XSDHayXDmBRR+TUNtX20g2OCYDO003iNERE0qvbdNjR5tMBLOqInNrvQ6+hEichIpKevsJkO+/V0ebTASzqiJxWzuUKJOfqoVbxFAkiIgA4eO4ShAA6tfdCsK+71HGajEUdkZOqW/U6IMof7b21EqchIpLevtqtTG5wwKFXgEUdkdP6MfkCAGBCPIdeiYiA3+fTOeLQK8Cijsgp5ZZUIim7BCoVMI5bmRARoaSiGqfz6+bT+UucpnlY1BE5oU21vXQDovwR5ON480aIiFrb/gzrfLrYIG+HbRdZ1BE5IdvQay8OvRIRAX8479VBe+kAFnVETievpBJHz1uHXnmKBBGRVd0iicGd2kucpPlY1BE5mZ/qhl47+iPIAZfsExG1tsLSKqQWlAJwzE2H67CoI3IydUXdrb3YS0dEBAC/phcBAOLCfOHvpZE4TfOxqCNyInkllThSN/TK+XRERACAvWnWoddhnR136BVgUUfkVH48/vuqV0fcLZ2IqLUJIbA3/SIAYFgsizoichA/HM8DAEzihsNERACAsxfLUGAwQuOqxoAox135CrCoI3Ia54srcCzHetbruJ4s6oiIAGBvmnU+3YCodnB3c5E4TcuwqCNyEj8kW3vpBscEINCHZ70SEQHA3tpFEsNiAyVO0nIs6oicxA/HrPPpJsaHSZyEiEgeTGYL9mdcAgDc6OCLJAAWdURO4ezFMpy8YICrWoVxcdzKhIgIAI5ll6DMWIN2nm7oEeordZwWY1FH5ATqeumGdW6Pdg68BxMRUWuqG3odEtsearVK4jQtx6KOSOGEEPi+dtUrh16JiH5Xt0jC0bcyqcOijkjhTl0oRXphGTSuaoyJC5Y6DhGRLJRWmXA0uwQAizoichAbj1l76W7uGgRfdzeJ0xARycOBjEswWwQ6Bngiwt9T6jitgkUdkYJZLALf1xZ1t/Xh0CsRUZ3ftzJRRi8dwKKOSNGOnL+M3JJKeGtdcXO3IKnjUK3du3dj0qRJCAsLg0qlwoYNG+p9XgiBefPmITQ0FB4eHhg9ejTS0tKkCUukUCzqiMihfJdk7aUbExfs8DulK0l5eTl69+6NxYsXN/j5d955B//5z3+wdOlSHDhwAF5eXhg7diyqqqrsnJRImfL1VUgvLINKBQyJUU5R5yp1ACJqGzVmC35Ktm5lcltvDr3Kyfjx4zF+/PgGPyeEwKJFi/Dyyy9j8uTJAIDPP/8cwcHB2LBhA6ZMmWLPqESKVNdLF99BB52ncuYas6eOSKF+PVuM4vJqBHhpMFRBwwtKl5mZifz8fIwePdr2nE6nw6BBg7Bv374GX2M0GmEwGOo9iOjqdp25CAC4sbPjHw32RyzqiBTqu6O5AIBbe4XCzYXf6o4iPz8fABAcXH/7meDgYNvn/mzBggXQ6XS2R0RERJvnJHJUZovAnjRrUTeiK4s6IpK5iuoabD5hLQBuT+DQq9LNnTsXer3e9sjOzpY6EpFsHcspQUmFCb7urugT4Sd1nFbFoo5IgbaeLEBFtRkR/h7oG9lO6jjUBCEh1rN5CwoK6j1fUFBg+9yfabVa+Pr61nsQUcN2pv4+9OqqsFEMZd0NEQEANtQOvd7RpwNUKsc/z9CZREdHIyQkBNu3b7c9ZzAYcODAAQwePFjCZETKUDefbrjChl4Brn4lUpyiMiN2155nODmhg8RpqCFlZWVIT0+3fZyZmYmkpCT4+/sjMjISzz33HN544w107twZ0dHReOWVVxAWFobbb79dutBEClBcZsTxnBIAwPAuLOqISOZ+OJYHs0UgPlyHmEBvqeNQAw4fPoyRI0faPp49ezYAYPr06VixYgX+/ve/o7y8HE888QRKSkowbNgwbN68Ge7u7lJFJlKEPWlFEALoHuqLYF/lfT+xqCNSmA21Gw7f3oe9dHI1YsQICCGu+nmVSoV//etf+Ne//mXHVETKZxt6VWAvHcA5dUSKkllUjqTsErioVZjEDYeJiGwsFoHdZ5S5lUkdFnVECrK+doHE0Nj2CPTRSpyGiEg+UvL0KC6vhrfWFf06KnNXABZ1RAphsQisO5IDALirL4deiYj+qG4rk6GxAYrdkF2Zd0XkhA5nXUbO5Up4a10xpkfD+5kRETmrXbah1yCJk7QdFnVEClHXS3drrxB4aFwkTkNEJB8lFdU4ev4yAOUukgBY1BEpQpXJjB+PXwAA3Nk3XOI0RETysietCBYBdAn2Rpifh9Rx2gyLOiIF2HqyAKXGGnTw88DAKH+p4xARyYrStzKpI9ui7rXXXoNKpar36Natm9SxiGSpbuj1zr4doFbzWDAiojoWi3CK+XSAzDcfjouLw7Zt22wfu7rKOi6RJApLq2zHgt3BY8GIiOpJydPjYqkRXhoX9I9S5lYmdWRdJbm6uiIkhKv4iK5l/ZFcmC0CfSP90InHghER1bPtVCEA4MbOgdC6KnsRmWyHXwEgLS0NYWFh6NSpE6ZOnYrz589LHYlIVoQQ+DbROvR6d78IidMQEcnP9lMFAIBR3ZU99ArIuKdu0KBBWLFiBbp27YoLFy5g/vz5uPHGG5GSkgIfH58GX2M0GmE0Gm0fGwwGe8UlksSxHD3SCsvg7qbGxN6hUschIpKVC/pKnMgzQKUCRnZjUSeZ8ePH2/4/Pj4egwYNQseOHfH111/j0UcfbfA1CxYswPz58+0VkUhy3yZmAwDGxYXA191N4jRERPLyy2nr0GtChB/aeyv/6ERZD7/+kZ+fH7p06YL09PSrXjN37lzo9XrbIzs7244JieyrymTGxqQ8ABx6JSJqyPba+XSjugdLnMQ+HKaoKysrw9mzZxEaevUhJq1WC19f33oPIqXaerIAhqoahOncMSQmQOo4RESyUlltxq/p1p0BRrOok9YLL7yAXbt24dy5c/jtt99wxx13wMXFBffff7/U0Yhk4ZvaBRJ39Qvn3nRERH+yN70IxhoLOvh5oEuwc+wMINs5dTk5Obj//vtRXFyMwMBADBs2DPv370dgoLJ3gyZqjNySSuxJs26meXc/HgtGRPRndateR3cPgkrlHL/4yraoW7t2rdQRiGTr28M5EAK4oZM/OgZ4SR2HiEhWLBZhWyThLPPpABkPvxJRwywWgW9qV73eN4ALJIiI/iwlT4/C2lMkBnVynvOwWdQROZjfzhYj53IlfNxdMb4n96YjIvqzulMkbuqi/FMk/ohFHZGD+eqwtZducp8wuLs5T2NFRNRYv58i4TxDrwCLOiKHUlJRjS0n8gEA9/WPlDgNEZH8/PEUiRFdnWtxJYs6Igey4Wguqmss6B7qi54duA8jEdGf1W047CynSPwRizoiByGEwNpDtQsk+oc7zRJ9IqKm2HrSOYdeARZ1RA4jKbsEp/NLoXVV444E7k1HRPRnhioTfjtrPUViXM8QidPYH4s6Igex9qC1l25Cr1DoPN0kTkNEJD87ThfCZBaIDfJGTKBznCLxRyzqiBxAaZUJG4/lAQCmDOQCCSKihtQtJBsX53y9dACLOiKHsPFYHipNZsQEemFAVDup4xARyU6VyYwdp63HJ45lUUdEcvXlwfMAgPsHRnKBBBFRA/akFaHSZEYHPw+n3R2ARR2RzKXk6pGSa4DGRY07+3KBBBFRQzanWIdex8QFO+0vvyzqiGRu9YEsAMDYniHw99JInIaISH5MZgu2n7ZuZeKsQ68AizoiWTNUmfBdknWBxIODuECCiKghBzMvoaTChAAvDQZE+UsdRzIs6ohk7LujuaioNiM2yBsDo523oSIiupa6Va+juwfDRe2cQ68Aizoi2RJCYPUB6wKJqYO4QIKIqCEWi/h9KxMn3HD4j1jUEcnUkfOXcTq/FO5uXCBBRHQ1x3JKUGAwwlvriiGxAVLHkRSLOiKZWr3f2kt3W+8w6Dx4ggQRUUM21/bSjegaCK2ri8RppMWijkiGLpVX44fkCwCAqYM6SpyGiEiehBDYksKh1zos6ohk6OvD2aiusaBnB1/Eh+ukjkNEJEupBaU4V1wBjYsaI7oGSR1HcizqiGTGbBFYtd+6N91DN0RxgQQR0VX8eNw6ojG8ayC8ta4Sp5Eeizoimdl1phA5lyuh83DDpN5hUschCZjNZrzyyiuIjo6Gh4cHYmJi8Prrr0MIIXU0ItkQQtiKuonxoRKnkQeWtUQy8/k+ay/dvf3D4aFx7km/zurtt9/GkiVLsHLlSsTFxeHw4cOYMWMGdDodnnnmGanjEcnCqQulyCgqh8ZVjVHdg6WOIwss6ohkJKu4HLvOXATABRLO7LfffsPkyZMxYcIEAEBUVBS+/PJLHDx4UOJkRPLxw3HraTsjOfRqw+FXIhlZfeA8hACGdwlEVHsvqeOQRIYMGYLt27fjzJkzAIBjx45h7969GD9+fIPXG41GGAyGeg8iJRNC4MfkuqFXTlOpw9KWSCYqq8346lA2AGDaDeylc2YvvfQSDAYDunXrBhcXF5jNZrz55puYOnVqg9cvWLAA8+fPt3NKIumcyDMgq7gC7m5q3NyNq17rsKeOSCY2JOVCX2lChL8HRrKRcmpff/01Vq9ejTVr1uDIkSNYuXIl3nvvPaxcubLB6+fOnQu9Xm97ZGdn2zkxkX19Xzv0enO3IHhx6NWGfxNEMiCEwIpfzwEApg+OcuoDqQl48cUX8dJLL2HKlCkAgF69eiErKwsLFizA9OnTr7heq9VCq9XaOyaRJOqveuXQ6x+xp45IBvZnXEJqQSk83FxwT/8IqeOQxCoqKqBW12+eXVxcYLFYJEpEJB/Hc/TIuVwJDzcXjOSGw/Wwp45IBlb8lgkAuLNvB57zSpg0aRLefPNNREZGIi4uDkePHsX777+PRx55ROpoRJKrWyAxqnsQt336ExZ1RBLLuVyBrScLAADTh0RJG4Zk4cMPP8Qrr7yCp59+GoWFhQgLC8OTTz6JefPmSR2NSFLccPjaWNQRSeyL/VmwCGBobAC6BPtIHYdkwMfHB4sWLcKiRYukjkIkK0ezS5BbUgkvjQvPem0A59QRSaiiugZfHjgPAHh4SLTEaaghLi4c3iGSi7peutE9guHuxu/NP2NRRySh/0vMgaGqBh0DPDGK25jIEs9bJZIHs0Xg+2PWrUwm9OLQa0OaXNRVVlYiNzf3iudPnDjRKoGInIXFIrC8dhuTGUOioOY2Jm3u5MmTeO+995CXZ/3BsH79+uu+RqXi+0IkB/vOFqOw1AidhxuHXq+iSUXdt99+i86dO2PChAmIj4/HgQMHbJ+bNm1aq4cjUrJdZy4io6gcPlpX3M1tTOxi/vz5uOWWWzB//nwcO3YMW7dulToSETXShiRrh9KE+FBoXDnQ2JAm/a288cYbSExMRFJSEpYvX45HH30Ua9asAcAhCqKm+uxX6zYmUwZG8DBqO/Hz80Pv3r2xdOlSLF++HMePH2/yn2EymbB161bs2bMHxcXFbZCSiP6symTG5pR8AMDtfTpInEa+mvSTxGQyITg4GADQr18/7N69G3fccQfS09M5REHUBGcKSrEnrQhqFfDQ4Cip4ziNW265BYB1SPXf//433nvvvSb/GXfeeSdCQ0Oxbt06tGvXDhUVFejVqxc2b97c2nGJqNa2UwUoM9agg58H+ndsJ3Uc2WpST11QUFC932z9/f2xdetWnDp1qlm/8RI5q//tyQAAjI0LQYS/p8RpnMfdd98NAJg5cyYuXryIF198scl/xvnz5/HJJ58gPDwcaWlp+Mc//oH4+PjWjkpEf7DhqHUe7OQ+YZx/fA1NKuq++OILBAXVn5yo0Wjw5ZdfYteuXa0ajEipCkurbA3UYzd2kjiNcxo/fjxuvfVWvPbaaygvL2/Sa93d3QFY277q6mrMmjULe/fubYuYRATgcnk1dqYWAgDuSODQ67U0qagLDw9HSEhIveeys7MBAEOHDm29VEQK9vlvWag2W9A30g/9OIwgicmTJ+PAgQMIDg7GkCFDsHTp0kafq/rMM8/g0qVLuOuuuzBz5kwsW7YMRUVFbZyYyHn9mHwBNRaBHqG+6MwN2q+pxctHunXrhnnz5qGioqI18hApWkV1DVYdyAIAPM5eOkm5uLhgwoQJeP755/Hyyy+jR48e+P777696/dmzZwEAU6dOhb+/P+bMmYObbroJp0+fxrfffmuv2EROZ8NR66pX9tJdX4uLuq1bt2LLli3o3LkzVqxY0QqRiJTr/xJzUFJhQqS/J8bEhVz/BdQmxo0bh8jISDzwwAM4fvw4PvzwQ6xevRobNmzAc8891+BrZs6ciejoaAwePBhPPvkkFi9ejJiYGLz88sucU0fURrIvVeBw1mWoVMCk3mFSx5G9Fu+jMGTIEBw4cACff/45/vnPf+LDDz/EokWLcOONN7ZGPiLFMFsElu21bmPy6LBouHCyr92dPXsWMTExWLhwIXr16nXFEWDLli1Dt27dGnxt3Z52b731Fg4dOoTc3Fxs3LgR27dvR1RUFNLT09s8P5Gz+a52b7ohMQEI0blLnEb+Wm1zrIceegh33303Fi5ciPHjx2PcuHF49913ER3N8yyJAODnE/k4V1wBnYcb7ukfLnUcpzRz5kykp6cjJCQE8fHx9R46nQ4A8NNPP13zz/j666+RlJRk+/jnn3/G6tWr2zI2kVMSQmB97dDrZO5N1yitviXzmDFj8Nhjj2H9+vXo0aMH/v73v6OsrKy1vwyRQxFCYOku65yshwZ3hKeGmw1LYevWrcjMzMSkSZNQWFiI3NxcvPHGG/D390dsbCwAoFOna891dHd3x8mTJ20fjxkzBikpKW2am8gZncgz4OzFcmhd1RjXk9NVGqPFP1mWLl2KQ4cO4dChQzh16hTUajV69uyJmTNnonfv3li7di169OiBdevWoX///q2RmcjhHMi8hGM5emhc1Zg+JErqOE6vJb1ty5Ytw3333YcRI0agT58+SE5O5ubrRG3g28QcAMDoHsHwdXeTOI1jaHFP3Ztvvgm9Xo+HHnoIO3bsQElJCRITE7F48WI88cQT+OWXXzBz5kw8/PDDLfo6CxcuhEqluuokZiI5+7i2l+6efuFo762VOA21pLctLi4OiYmJuPHGG3Hu3Dl07NgRmzZtaquoRE7JWGO2nfV6Tz9OV2msFvfU1e1Tdy2PPvooXnnllWZ/jUOHDuHjjz/mCjNySKfzDdiRehFqFbcxkYuW9rZpNBrce++9uPfee9swJZHz2n6qECUVJoT4uuPGzoFSx3EYrT6nriFBQUH45ZdfmvXasrIyTJ06FZ9++inateNGreR4PtltPRJsfM9QRLX3kjgNAU3rbRNC2DkdEX1z2NphdGffDtwpoAnsMltbpVJh+PDhzXrtrFmzMGHCBIwePRpvvPFGKycjals5lyuwMcl6JNgTN7GXTk4a29vW2JMmiKh1FBiqsOvMRQDA3Rx6bRJZL8Fbu3Ytjhw5gkOHDjXqeqPRCKPRaPvYYDC0VTSiRvl0dwZqLAJDYwPQO8JP6jhERLK37kguLALo37EdOgV6Sx3Hodhl+LU5srOz8eyzz2L16tW2A7SvZ8GCBdDpdLZHREREG6ckurqiMiPWHrIOITw9IlbiNERE8ieEsA29cj/PppNtUZeYmIjCwkL07dsXrq6ucHV1xa5du/Cf//wHrq6uMJvNV7xm7ty50Ov1tkdjFnEQtZXlv2bCWGNB73AdhsQESB2HiEj2jpy/jIyicni4uWBCPI8FayrZDr+OGjUKycnJ9Z6bMWMGunXrhjlz5lxxvA8AaLVaaLXcLoKkV1plwuf7sgAAT42I5T5mRESN8M1h695043uFwFsr2xJFtmT7N+bj44OePXvWe87LywsBAQFXPE8kN6v2n0dpVQ1ig7wxpkew1HGIiGSvoroGPxy/AAC4px+nTzWHbIdfiRxVZbUZ/9tj3cZk5vAYqLkcn4joujan5KPMWIMIfw8MivaXOo5Dkm1PXUN27twpdQSi61pz8DyKy6sR4e+ByX04J4SIqDHqhl7v7hvBX4abiT11RK2oymTGJ7utR4I9NTwWbi78FiMiup5zReXYl1EMlQq4q18HqeM4LP7EIWpF3ybmoMBgRKjOnQ0TEVEjfXnwPABgeJdAhLfzlDiN42JRR9RKTGYLluy09tI9eVMnaF2vXKFNRET1GWvM+CbROvT6wMBIidM4NhZ1RK1k/ZFc5JZUor23FlPYMBERNcqWEwW4VF6NEF933NwtSOo4Do1FHVErMJkt+HBHGgDgiZui4e7GXjoiosZYc8C6p+e9AyLgynnILcK/PaJWsP5oLrIvVSLAS4MHb+godRwiIodw9mIZ9mdcgloFTBnAvelaikUdUQvVmC1YvCMdAPDk8E7w1DjUTkFERJJZW7tAYmTXIIT5eUicxvGxqCNqofVHc5FVXMFeOiKiJqgymfFt3QKJQZyH3BpY1BG1wB976Z64ib10RESNteVEPi5XmBCmc8eIrlwg0RpY1BG1wLqjuThXXAF/9tIRETXJ6gPWodf7BkTChSdItAoWdUTNVF1jwX+2W1e8zhzeCV5a9tIRETVGemEpDmZegotahfu4QKLVsKgjaqZvE3OQc9m6L920G6KkjkNE5DBW7bf20t3cLQghOneJ0ygHizqiZjDWmPHRL9ZeuqdHxMBDw33piIgao8xYY1sg8dBgTltpTSzqiJph7cFs5OmrEOLrzlVb1CZyc3Px4IMPIiAgAB4eHujVqxcOHz4sdSyiFvu/xByUGWsQE+iFYbHtpY6jKJwERNREldVm24rXWTfH8vQIanWXL1/G0KFDMXLkSGzatAmBgYFIS0tDu3btpI5G1CIWi8DKfecAANOHREGl4gKJ1sSijqiJPt93DoWlRnTw88C9/cOljkMK9PbbbyMiIgLLly+3PRcdHS1hIqLWsTe9CBkXy+GtdcWdfdl+tjYOvxI1gaHKhCW7zgIAnr+lC7Su7KWj1rdx40b0798f99xzD4KCgpCQkIBPP/30qtcbjUYYDIZ6DyI5WvnbOQDA3f3C4c0dA1odizqiJvjfnkyUVJgQG+SNOxI6SB2HFCojIwNLlixB586dsWXLFjz11FN45plnsHLlygavX7BgAXQ6ne0REcEtIkh+sorL8UtqIQAukGgrLOqIGqm4zIhlezIAAH+7pQs3y6Q2Y7FY0LdvX7z11ltISEjAE088gccffxxLly5t8Pq5c+dCr9fbHtnZ2XZOTHR9n+/LghDA8C6B6BToLXUcRWJRR9RIi3ecRXm1Gb066DCuZ4jUcUjBQkND0aNHj3rPde/eHefPn2/weq1WC19f33oPIjkpN9bg68PWXzYeHhIlbRgFY1FH1Ag5lyuwan8WAODFsV25Yova1NChQ5GamlrvuTNnzqBjRw5ZkWNafzQXpVU1iArwxPAugVLHUSwWdUSN8P7WM6g2WzAkJgA3dua+StS2nn/+eezfvx9vvfUW0tPTsWbNGnzyySeYNWuW1NGImkwIYVsgMW1wFNScutJmWNQRXcepCwasP5oLAJgzrht76ajNDRgwAOvXr8eXX36Jnj174vXXX8eiRYswdepUqaMRNdnutCKkFZbBU+OCe7gNVJviemKi63h3SyqEACbEh6J3hJ/UcchJTJw4ERMnTpQ6BlGLfbrbusBsyoBI+Lq7SZxG2dhTR3QNBzKK8cvpQriqVXhhTFep4xAROZQTeXrsTS+Ci1qFGUOjpI6jeCzqiK7CYhF466dTAID7BkQgur2XxImIiBzL//ZkAgBu7RWKCH9PidMoH4s6oqv4IfkCjuXo4aVxwXOju0gdh4jIoeSVVOL7Y3kAgMdv5DF39sCijqgBxhoz3tl8GgDw5PAYBPpoJU5ERORYVvx2DjUWgRs6+SM+3E/qOE6BRR1RA77Yl4Wcy5UI8tHiMf6GSUTUJKVVJnx5wLpZ9hM3dZI4jfNgUUf0JyUV1fjwl3QAwAtjusJTw0XiRERN8dWhbJQaaxAb5I0RXYKkjuM0WNQR/cl/tqdDX2lC12Af3NWPeyoRETWFyWzBZ3utCyQevzGamw3bEYs6oj/IuFiGz/edAwC8PLE7XNgYERE1yU/JF5Cnr0J7by0m9+kgdRynwqKO6A8WbDqNGovAzd2CcGNnnk9IRNQUFovAkp1nAQDTB3eEu5uLxImcC4s6olq/nS3C1pMFcFGr8I9bu0kdh4jI4Ww/XYjT+aXw1rriocFRUsdxOizqiACYLQKv/2DdaPjBQZGIDfKROBERkWMRQuCjHdZFZg/e0BE6Tx4JZm8s6ohgXal16oIBvu6ueJYbDRMRNdmv6cU4ll0Crasajw7jVlBSYFFHTk9fYcJ7P6cCAJ6/pQv8vTQSJyIicjwf7UgDANw/MJIbtkuERR05vQ+2p+FSeTU6B3njwRs6Sh2HiMjhJGZdwv6MS3BzUXGzYQmxqCOnll5YatvCZN6kHnBz4bcEEVFTfVS7YfudCeEI8/OQOI3z4k8wclpCCMz//iRqLAKjuwdzCxMiomZIydVjR+pFqFXAUyNipI7j1FjUkdPacqIAe9KKoHFR4+UJ3aWOQ0TkkP6709pLNzE+DFHtvSRO49xY1JFTqqw24/UfTgIAnhzeiQ0REVEzpBeWYlNKPgBg1shYidMQizpySkt2piO3pBId/Dzw9Ag2REREzfHvbWkQAhjTIxhdQ7i/p9RY1JHTySoux9LdGQCAlyd0h4eGx9gQETXV6XwDfjx+AYB1OyiSHos6cipCCLy28QSqaywYFtse43qGSB2JiMgh/XvrGQDAhF6h6B7qK3EaAljUkZPZcqIAO1Ivws1FhfmT46BSqaSORETkcFJy9dhyogAqFfDc6M5Sx6FaLOrIaZQba/Cv708AAJ68KQYxgd4SJyIickx1vXSTe4ehczDn0smFbIu6JUuWID4+Hr6+vvD19cXgwYOxadMmqWORA/vPL2nI01chvJ0HV2kRETXT0fOXsf10IVzUKp6VLTOyLerCw8OxcOFCJCYm4vDhw7j55psxefJknDhxQupo5IBS80uxbE8mAGD+bXFcHEFE1Ezv1/bS3ZHQAdHcDkpWXKUOcDWTJk2q9/Gbb76JJUuWYP/+/YiLi5MoFTkii0XgH+uTUWMRuKVHMEZ1D5Y6EhGRQzp07hL2pBXBVa3Cs6M4l05uZFvU/ZHZbMY333yD8vJyDB48+KrXGY1GGI1G28cGg8Ee8Ujmvjx0HolZl+GlccH82/gLARFRcwgh8N6WVADAPf0jEOHvKXEi+jPZDr8CQHJyMry9vaHVajFz5kysX78ePXr0uOr1CxYsgE6nsz0iIiLsmJbkqLC0Cgs3nQYA/G1MVx40TUTUTDtTL+JA5iVoXNX4682clyxHsi7qunbtiqSkJBw4cABPPfUUpk+fjpMnT171+rlz50Kv19se2dnZdkxLcvT6D6dQWlWDXh10mD4kSuo4REQOyWwRtl+QZwyJ4i/IMiXr4VeNRoPYWOtvA/369cOhQ4fwwQcf4OOPP27weq1WC61Wa8+IJGO/nC7A98fyoFYBb93RCy5q7klHRNQc647kILWgFL7urnhqRIzUcegqZN1T92cWi6XenDmiqykz1uDl9SkAgEeHRaNXuE7iREREjqnKZLateJ01MhZ+nhqJE9HVyLanbu7cuRg/fjwiIyNRWlqKNWvWYOfOndiyZYvU0cgBvLv5NPL0VYj098TsW7pKHYeIyGGt/O0cLuirEKZz5zQWmZNtUVdYWIiHHnoIFy5cgE6nQ3x8PLZs2YJbbrlF6mgkc4lZl/H5/iwA1mFX7klHRNQ8JRXVWLwjHQAwe0xXuLuxPZUz2RZ1y5YtkzoCOaAqkxlz/u84hADu6ReOYZ3bSx2JiMhh/XfnWRiqatAtxAd3JHSQOg5dh0PNqSO6nv9sT0N6YRkCfbT454TuUschInJYuSWVWPHbOQDAnHHduNjMAbCoI8VIztHj490ZAIA3bu/JybxERC3w9qbTqK6x4IZO/hjRNVDqONQILOpIEaprLHjx22MwWwQmxodibFyI1JGIiBzW4XOXsPFYHlQq4OUJPaBSsZfOEbCoI0X4aEc6TueXwt9Lw6PAiIhawGIRmP+9daP/KQMi0LMDt4RyFCzqyOGl5Optq7P+NTkOAd7cgJqIqLm+PZKD5Fw9fLSu+NsYbgnlSFjUkUMz1pgx++skmC0CE+JDMTE+TOpIRK1u4cKFUKlUeO6556SOQgpXWmXCO5tTAQDPjOqM9vwl2aGwqCOH9sG2NJwpKEN7bw1en9xT6jhEre7QoUP4+OOPER8fL3UUcgKLd5xFUZkR0e29uNGwA2JRRw7ryPnLWLrrLADgjdt7wd+Lq11JWcrKyjB16lR8+umnaNeundRxSOGyisvx2d5MAMDLE7pD48oSwdHwHSOHVFFdg9lfJcEigNv7hGFcT652JeWZNWsWJkyYgNGjR1/zOqPRCIPBUO9B1FRv/ngK1WYLbuzcHjd3C5I6DjWDbE+UILqWN388hXPFFQjVuWM+h11JgdauXYsjR47g0KFD1712wYIFmD9/vh1SkVLtOF2In08WwEWtwryJ3MLEUbGnjhzOjtRCrD5wHgDw3j29ofNwkzgRUevKzs7Gs88+i9WrV8Pd3f2618+dOxd6vd72yM7OtkNKUorKajPmbUwBADw6LBqdg30kTkTNxZ46ciiXyqvx92+PAwBmDI3C0Fie7UrKk5iYiMLCQvTt29f2nNlsxu7du/HRRx/BaDTCxeX3g9W1Wi20Wq5SpOZZvCMd2ZcqEapzx7OjOksdh1qARR05DCEE5vzfcVwsNSI2yBtzxnWTOhJRmxg1ahSSk5PrPTdjxgx069YNc+bMqVfQEbXE2Ytl+Hi3dcHZq5N6wEvLssCR8d0jh/HlwWxsPVkANxcVPpjSB+5u/MFGyuTj44OePevPFfXy8kJAQMAVzxM1lxACr2xIgcksMLJrII9XVADOqSOHcPZiGV7/wXpszd/HdkNcGI+tISJqiY3H8vDb2WJoXdWYf1tPLo5QAPbUkewZa8x4du1RVJrMGBobgEeHRUsdicjudu7cKXUEUhB9pQmv/3AKAPDXm2MRGeApcSJqDeypI9l7Z3MqUnINaOfphv93Tx+o1fxtkoioJd7ZfBpFZUZ0CvTC4zd1kjoOtRIWdSRrO04XYlntDufv3t0bIbrrb+9ARERXt+9ssW1bqDdu7wmtK+cnKwWLOpKtQkMVXvjmGADg4SFRGN0jWOJERESOrbLajJfWWbeFun9gJIbEcFsoJWFRR7Jktgg8uzYJxeXV6B7qi5fGc/sSIqKW+ve2M8gqrkCIrzvm3sp2VWlY1JEsffhLGvZlFMNT44KPHkjg9iVERC10LLsE/9uTAQB4686e8HXnaTxKw6KOZOe39CJ8sD0NAPDmHT0RE+gtcSIiIsdWXWPB3789DosAJvcJw83dOJ1FiVjUkawUllbh2a+SIARwb/9w3JEQLnUkIiKHt3hHOlILShHgpcGrk+KkjkNthEUdyUaN2YJnvjyKi6VGdAn2xvzbuHM+EVFLncjT47870wEAr90WB38vjcSJqK2wqCPZ+Pe2M9ifcQmeGhf8d2o/eGg4j46IqCWqTGY8/1USTGaBsXHBmBgfKnUkakMs6kgWfjldgMU7rIdKv31XPGKDOI+OiKil3tuSijMFZWjvrcFbd/TiUWAKx6KOJHe+uALPf/X7fnSTeodJnIiIyPH9drYIy361bt7+9l3xCPDWSpyI2hqLOpJUZbUZT65KhL7ShD4RfvjHrd2ljkRE5PAMVSa88PUxCAHcPzACo7pztaszYFFHkhFC4J8bknHqggHtvTVY8mBfaFz5T5KIqKVe23gCefoqRPp74uUJPaSOQ3bCn6AkmVX7s7DuSC5c1Cp8eH9fhOo8pI5EROTwNiVfwLojuVCrgPfv7Q0vravUkchOWNSRJA5kFGP+9ycBAHPGdcXgmACJExEROb6cyxWY83/Ws11nDo9B/yh/iRORPbGoI7vLK6nE06uPoMYicFvvMDx+YyepIxEROTxT7V6fhqoa9I7ww3Oju0gdieyMRR3ZVZXJjCe+OIzi8mr0CPXF23fFc4k9EVEr+PfWMzhyvgQ+Wld8OCWBc5SdEN9xshshBOb833Gk5Brg76XBJw9xg2EiotawJ+0iluyy7vW58K54RAZ4SpyIpMCijuzmvzvP4rukPLiqVVj8QF+Et2OjQ0TUUoWlVXi+9szsBwZFYgJPjXBaLOrILraeLMB7P6cCsJ49yIURREQtZ7EIzP7qGIrKqtEtxAfzJnL7EmfGoo7a3KkLBjy39iiEAKbd0BEP3tBR6khERIrwwfY07E0vgoebCz56IAHubpzS4sxY1FGbKiytwqMrDqG82owhMQGYN4m/RRIRtYZfThfgg+1pAIA3bu+J2CAfiROR1FjUUZupMpnx+OeJyNNXoVOgF5ZM7Qc3F/6TIyJqqXNF5XhubRIA6wjIXf3CpQ1EssCfsNQmLBaBv319DMeyS+Dn6YbPpg+AztNN6lhERA6voroGM1clwlBVg76RfniF8+ioFos6ahNvbzmNH5MvwM1FhY8f7Ieo9l5SRyIicnhCCMxdl4zT+aVo763Bf6f24350ZMN/CdTqVu3Pwse7MgAA79wdj0GduNKViKg1rPjtHL5LyoOLWoWPHuiLEJ271JFIRljUUav65XQB5n2XAgCYfUsX3JHAeR5ERK1hT9pFvPHjKQDA3PHdcAN/YaY/YVFHreZYdglmrT4KiwDu7heOv94cK3UkIiJFSC8sw9Orj8BsEbgzoQMeHRYtdSSSIRZ11CrOFZXjkRWHUGky48bO7fHWHb14pisRUSu4XF6NR1ceQmlVDfp3bIcFd7F9pYaxqKMWKyozYvrygygur0ZcmC+WPMiJu0REraG6xoKnViciq7gCHfw8sHRaP2hducEwNYw/ealFSqtMeHj5QWQVVyC8nQeWzxgAb62r1LGIiByeEAKvbkzB/oxL8NK4YNnD/dHeWyt1LJIx2RZ1CxYswIABA+Dj44OgoCDcfvvtSE1NlToW/UGVyYwnv0hESq4B/l4afP7IQAT5cCUWEVFr+HRPBr48mA2VCvjP/QnoFuIrdSSSOdkWdbt27cKsWbOwf/9+bN26FSaTCWPGjEF5ebnU0QiA2SIw++sk/Ha2GF4aF6ycMRCdAr2ljkVEpAjfJeXirZ9OAwD+eWt3jOoeLHEicgSyHSfbvHlzvY9XrFiBoKAgJCYm4qabbpIoFQHWIYF/rEvGT8n50Lio8elD/dErXCd1LCIiRfg1vQgvfHMMAPDI0GiudKVGk21R92d6vR4A4O/vf9VrjEYjjEaj7WODwdDmuZyNEAJv/ngKXx3OhloFLJrSB0Ni20sdi4hIEU7k6fHkF4kwmQUmxIfi5QndudKVGk22w69/ZLFY8Nxzz2Ho0KHo2bPnVa9bsGABdDqd7REREWHHlM7ho1/S8b+9mQCAhXfG49ZeoRInIiJShuxLFXh4+SGUGWtwQyd/vH9vb6jVLOio8RyiqJs1axZSUlKwdu3aa143d+5c6PV62yM7O9tOCZ3D//Zk4P9tPQMAeGViD9w7gEUzEVFrqNsa6mKpEd1CfPDxtP7cuoSaTPbDr3/5y1/www8/YPfu3QgPv/aRU1qtFlotl3u3hVX7s2zH0zw3ujPneBARtRJ9hQnTlh1ExsVyhOncsXzGAOg83KSORQ5ItkWdEAJ//etfsX79euzcuRPR0SwipPJtYg5e3mA9z3Xm8Bg8O6qzxImIiJShzFiDh1ccxKkLBrT31mLVY4MQqvOQOhY5KNkOv86aNQurVq3CmjVr4OPjg/z8fOTn56OyslLqaE5lw9FcvPitdRXWw0OiMGdcV07aJWpj3KfTOVSZzHhs5SEcPV8CnYcbVj3GraGoZWRb1C1ZsgR6vR4jRoxAaGio7fHVV19JHc1pbDyWh9lfJ0EI4P6BkZg3sQcLOiI74D6dylddY8FTqxKxP+MSvLWu+PyRgdxcmFpM1sOvJJ0fjufh+a+SYBHAlAERePP2nlyFRWQn3KdT2aprLPjLmiPYkXoR7m5qLJveH70j/KSORQog26KOpLPxmLWgM1sE7ukXjrfu6MWCjkhC19unk3t0Og5jjRmzVh/BtlOF0Liq8fG0/hjUKUDqWKQQsh1+JWlsOJqL59YehdkicG//cLx9VzwLOiIJNWafTu7R6RiqTGbM/CIR204VQutqPY1neJdAqWORgrCoI5tvDmdj9tfWIdf7B0Zg4Z0s6Iik1ph9OrlHp/xVmcx48ovEPwy5DmBBR62Ow68EAPhifxZeqd22ZOqgSLw+mXPoiKTW2H06uUenvFVU1+DJLxKxJ60IHm4uWPZwfwyJ4fGK1PpY1BH+tyfDtrHwjKFRXOVKJDHu06kc+goTHll5CIlZl+GpccFnDw/ADZxDR22ERZ0TE0Jg0bY0fLA9DQDw9IgYvDiW+9ARSW3WrFlYs2YNvvvuO9s+nQCg0+ng4cGNaR1FoaEK05YdRGpBKXzdXbF8xgD069jwYhei1sCizkkJIfD6D6fw2a+ZAIC/3dIFf7k5lgUdkQwsWbIEADBixIh6zy9fvhwPP/yw/QNRk2UVl+PBZQeQfakSQT5afP4o96GjtseizgnVmC2Yuy4Z3yTmAABem9QDDw/l8A6RXHCfTsd2Ms+A6csP4mKpER0DPPHFI4MQGeApdSxyAizqnEyVyYy/rDmKbacKoFYB79zdG3f3u/oEbCIiarzdZy7i6dVHUGasQfdQX6x8ZACCfNyljkVOgkWdEzFUmfDYysM4mHkJGlc1Fj/QF7f0CJY6FhGRInx58Dxe3pACs0VgULQ/PnmoP3QeblLHIifCos5J5Our8PDygzidXwofrSv+N527mBMRtQaLReDdn1OxZOdZAMCdCR2w8K54aFy5FSzZF4s6J5BeWIqHlh1Enr4KgT5arJgxAHFhOqljERE5vCqTGS98cww/HL8AAHh2VGc8N7ozF52RJFjUKdzBzEt4/PPD0Fea0CnQCytnDESEPyfsEhG1VF5JJZ78IhHJuXq4uaiw4M54zlEmSbGoU7CNx/LwwtfHUG22ICHSD8umD4C/l0bqWEREDu/QuUt4alUiisqq0c7TDYun9uUpESQ5FnUKJITAkl1n8c7mVADA2LhgfDAlAe5uLhInIyJyfGsOnMerG1NgMgt0C/HBpw/15wgIyQKLOoWprrHg5Q3J+PqwdQ+6R4ZG458TusOF57gSEbVIlcmMf/1wEmsOnAcATOgVinfviYenhj9KSR74L1FBSiqqMXNVIvZnXIJaBcybyE2FiYhaQ1ZxOZ5efQQn8gxQqYAXxnTF0yNiuCCCZIVFnUKkF5bh8c8PI7OoHF4aF3z0QF+M7BYkdSwiIoe3OeUCXvzmOEqNNWjn6YZ/39cHI7qyfSX5YVGnALvOXMRf1hxBaVUNOvh54H/T+6N7KM8YJCJqieoaCxZsOoXlv54DAPTr2A4f3p+AMD8PaYMRXQWLOgcmhMCyvZl466dTsAigf8d2WDqtH9p7a6WORkTk0NILy/DcV0eRkmsAADxxUye8OLYr3Fy4oTDJF4s6B1VlMuMf65Kx7mguAODufuF4846e0LpyhSsRUXMJIbDqwHm8+eNJVJks8PN0wzt3xWNMXIjU0Yiui0WdA8orqcTMVYk4nqOHi1qFf97aHTOGRnHCLhFRCxSVGTHn2+PYfroQAHBj5/Z4757eCPZ1lzgZUeOwqHMwv6UX4S9fHsWl8mr4ebrhvw/0xZBYbnhJRNQSm5Iv4JXvUlBUVg2NixpzxnfDjCFRUHM7KHIgLOochBACn+zOwNubT8MigLgwXyx9sB83vCQiaoGiMiPmfZeCn5LzAQBdg32waEofLjYjh8SizgHoK0148Ztj+PlkAQDgrr7W+XM8IYKIqHmEENh4LA+vbTyByxUmuKhVeHpEDP5ycyznJpPDYlEncyfy9Hh69RFkFVdA46LGK5N64MFBkZw/R0TUTDmXK/Dqdydsc+e6h/ri3bvj0bODTuJkRC3Dok6mhBBYc/A85n9/EtU1FnTw88B/p/ZF7wg/qaMRETkkk9mCZXsz8cG2NFSazHBzUWHWyFg8PSIWGlduVUKOj0WdDJVWmfCP9Sn4/lgeAODmbkF4/97e8PPUSJyMiMgxHTp3Cf9cn4wzBWUAgIHR/njz9p7oHOwjcTKi1sOiTmaOZZfgmbVHkVVcAVe1Cn8f1xWPDevEFVhERM1wQV+JhZtO47sk6y/J/l4a/OPW7rirbwdOYyHFYVEnExaLwP/2ZuCdzamosQh08PPAf+5PQL+O7aSORkTkcCqrzfhkdwaW7EpHlckClQq4r38E5ozrhnZeHPUgZWJRJwMFhir87etj2JteBAC4tVcIFtwRD52nm8TJiIgci8Ui8P3xPLy96TTy9FUArEcovjopDr3CuRCClI1FncS2nMjHnP87jpIKE9zd1Jg3MQ73D4zgsAARURMIIbA7rQhvbzqNkxes57WG6dwx99bumBgfyjaVnAKLOomUVpnw+g8n8fXhHABAzw6+WHRfAmKDvCVORkTkWJKyS/D2ptPYl1EMAPDWuuLJmzrhsRs7wUPDPefIebCok8DBzEuY/XUSci5XQqUCnripE/52S1cuqSciaoLkHD0+2H4G205Z95vTuKgxbXBHzBoZC3/OmyMnxKLOjqpMZry7JRWf/ZoJIYDwdh54/94+GBjtL3U0IiKHcSy7BB9sT8MvtZsHq1XAHQnheP6Wzghvx6MTyXmxqLOTI+cv44VvjiHjYjkA4N7+4XhlYg/4uHMxBBFRYyRmXcZHv6RhR+pFANZi7vY+HTDr5ljEBHLqChGLujZWWW3G//s5Fctqe+eCfLRYeFcv3NwtWOpoRESyJ4TAztSLWLLrLA5mXgJQW8wldMBfRsaiE4s5IhsWdW1of0Yx5q5LRmaRtXfuzoQOmDepB0+GICK6jiqTGd8l5eJ/ezKRVmg9BcLNRYXb+3TA0yNjEd3eS+KERPLDoq4N6CtNWLjpNL48eB4AEOyrxYI72TtHRHQ9+foqrNqfhTUHz+NSeTUA62rW+wdG4JFh0QjVeUickEi+WNS1IiEENqXk47WNJ1BYagQA3D8wEi+N7wadB+fOERE1RAiBfRnFWLU/Cz+fKECNRQCw7jM3Y2g07hsYAV/OPya6LhZ1rSS3pBLzNqRge+1qrOj2XlhwZy/c0ClA4mRERPJUXGbEuiO5WHvoPM7WLiIDgIHR/pgxJAq39AiGqwu3eiJqLBZ1LWQyW7BsbyY+2JaGSpMZbi4qPDU8Bk+PjIW7Gze9JCL6oxqzBXvSivBtYg5+PpkPk9naK+elccHtCR3w4A0d0T3UV+KURI6JRV0L7M8oxisbUmyTeAdG+ePNO3qic7CPxMmIiOTlZJ4B64/m4LukPNv0FACID9dhyoBI3NYnDN5a/kgiagl+BzVDvr4Kb/10ChuP5QEAArw0mHtrd9zVtwPPFyQiqpVVXI4fjl/AxqQ8pBaU2p7399Jgcp8w3NMvAj3C2CtH1FpY1DWBscaMZXszsfiXdJRXm6FSAVMHReKFMV25TQkREYDzxRX4MfkCNqVcwPEcve15jYsao7oH4faEDhjZNYjHIhK1ARZ1jSCEwLZThXjjx5PIKq4AAPSN9MO/JvdEzw46idMREUlHCIFTF0rx88l8bE7Jx+n833vk1CpgSEx73NY7DGN7hnAXAKI2xqLuOk7mGfDmTyfxa3oxAOuJEC+N74bb+3SAWs2hViJyPlUmM/ZlFGPH6UJsO1mAPH2V7XNqFXBDpwDc2isUY+NCEOijlTApkXORdVG3e/duvPvuu0hMTMSFCxewfv163H777Xb52oWGKry/9Qy+OpwNIQCNqxqPDYvGrJGx8OJkXiKyg8WLF+Pdd99Ffn4+evfujQ8//BADBw60ew4hBDKKyrHnzEXsOnMRv50thrHGYvu8u5saw2LbY0xcCEZ3D4a/F6ejEElB1tVJeXk5evfujUceeQR33nmnfb6msQaf7M7AJ7szUGkyAwAmxodizrhuiPD3tEsGIqKvvvoKs2fPxtKlSzFo0CAsWrQIY8eORWpqKoKCgtr86xcYqrDvbDF+TS/Cb2eLkVtSWe/zoTp3jOgahNHdgzAkpj08NNzCiUhqsi7qxo8fj/Hjx9vla1XXWLD20Hn8Z3saisqsR9P0jfTDPyd0R7+O/nbJQERU5/3338fjjz+OGTNmAACWLl2KH3/8EZ999hleeumlVv96uSWVOJhZjIOZl3EgoxgZReX1Pq9xUaN/VDvc1CUQI7sGoUuwN1f7E8mMrIu6pjIajTAaf9//yGAwNOp1204W4PU/LIKICvDEnHHdMK5nCBstIrK76upqJCYmYu7cubbn1Go1Ro8ejX379l1xfXPbPgCYuy4Zu1IL682LAwCVCugZpsOQmAAMiW2PgVH+7I0jkjlFFXULFizA/Pnzm/y6fEMVsoor0N5bi2dHd8aUARFw49E0RCSRoqIimM1mBAcH13s+ODgYp0+fvuL65rZ9gLWHLk9fBRe1Cj3DfDEgyh83dArAgCh/6Dy5WpXIkSiqqJs7dy5mz55t+9hgMCAiIuK6r7tvQASMNRZMGRDBRRBE5HCa2/YBwF9GxmLmTZ3QJ9IPnhq2f0SOTFHfwVqtFlpt05fPu7mo8eiw6DZIRETUdO3bt4eLiwsKCgrqPV9QUICQkJArrm9u2wcAA6M5Z5hIKTjGSEQkMxqNBv369cP27dttz1ksFmzfvh2DBw+WMBkRyZmse+rKysqQnp5u+zgzMxNJSUnw9/dHZGSkhMmIiNrW7NmzMX36dPTv3x8DBw7EokWLUF5eblsNS0T0Z7Iu6g4fPoyRI0faPq6bMzJ9+nSsWLFColRERG3vvvvuw8WLFzFv3jzk5+ejT58+2Lx58xWLJ4iI6qiEEELqEG3FYDBAp9NBr9fD19dX6jhE1Ej83m0Z/v0ROa6WfP9yTh0RERGRArCoIyIiIlIAFnVERERECsCijoiIiEgBWNQRERERKQCLOiIiIiIFYFFHREREpAAs6oiIiIgUgEUdERERkQKwqCMiIiJSABZ1RERERArgKnWAtlR3rK3BYJA4CRE1Rd33rIKPpm5TbPuIHFdL2j9FF3WlpaUAgIiICImTEFFzlJaWQqfTSR3D4bDtI3J8zWn/VELBvwpbLBbk5eXBx8cHKpXqmtcaDAZEREQgOzsbvr6+dkrY9nhfjkep99aU+xJCoLS0FGFhYVCrOUukqdj2WSn13pR6X4By781e7Z+ie+rUajXCw8Ob9BpfX19F/UOqw/tyPEq9t8beF3vomo9tX31KvTel3heg3Htr6/aPvwITERERKQCLOiIiIiIFYFFXS6vV4tVXX4VWq5U6SqvifTkepd6bUu/L0Sn5fVHqvSn1vgDl3pu97kvRCyWIiIiInAV76oiIiIgUgEUdERERkQKwqCMiIiJSAMUWdYsXL0ZUVBTc3d0xaNAgHDx48JrXf/PNN+jWrRvc3d3Rq1cv/PTTT/U+L4TAvHnzEBoaCg8PD4wePRppaWlteQtX1ZR7+/TTT3HjjTeiXbt2aNeuHUaPHn3F9Q8//DBUKlW9x7hx49r6Nq7QlPtasWLFFZnd3d3rXSOX96wp9zVixIgr7kulUmHChAm2a+Tyfu3evRuTJk1CWFgYVCoVNmzYcN3X7Ny5E3379oVWq0VsbCxWrFhxxTVN/d6lKym1/VNq2wew/QMcp/2TddsnFGjt2rVCo9GIzz77TJw4cUI8/vjjws/PTxQUFDR4/a+//ipcXFzEO++8I06ePClefvll4ebmJpKTk23XLFy4UOh0OrFhwwZx7Ngxcdttt4no6GhRWVlpr9sSQjT93h544AGxePFicfToUXHq1Cnx8MMPC51OJ3JycmzXTJ8+XYwbN05cuHDB9rh06ZK9bkkI0fT7Wr58ufD19a2XOT8/v941cnjPmnpfxcXF9e4pJSVFuLi4iOXLl9uukcP7JYQQP/30k/jnP/8p1q1bJwCI9evXX/P6jIwM4enpKWbPni1OnjwpPvzwQ+Hi4iI2b95su6apf190JaW2f0pt+4Rg+1fHUdo/Obd9iizqBg4cKGbNmmX72Gw2i7CwMLFgwYIGr7/33nvFhAkT6j03aNAg8eSTTwohhLBYLCIkJES8++67ts+XlJQIrVYrvvzyyza4g6tr6r39WU1NjfDx8RErV660PTd9+nQxefLk1o7aJE29r+XLlwudTnfVP08u71lL369///vfwsfHR5SVldmek8P79WeNadj+/ve/i7i4uHrP3XfffWLs2LG2j1v690XKbf+U2vYJwfbvahyh/ZNb26e44dfq6mokJiZi9OjRtufUajVGjx6Nffv2Nfiaffv21bseAMaOHWu7PjMzE/n5+fWu0el0GDRo0FX/zLbQnHv7s4qKCphMJvj7+9d7fufOnQgKCkLXrl3x1FNPobi4uFWzX0tz76usrAwdO3ZEREQEJk+ejBMnTtg+J4f3rDXer2XLlmHKlCnw8vKq97yU71dzXe/7rDX+vpydUts/pbZ9ANu/a1FK+2fPtk9xRV1RURHMZjOCg4PrPR8cHIz8/PwGX5Ofn3/N6+v+25Q/sy00597+bM6cOQgLC6v3j2fcuHH4/PPPsX37drz99tvYtWsXxo8fD7PZ3Kr5r6Y599W1a1d89tln+O6777Bq1SpYLBYMGTIEOTk5AOTxnrX0/Tp48CBSUlLw2GOP1Xte6verua72fWYwGFBZWdkq/76dnVLbP6W2fQDbv6tRUvtnz7bPtcVpyWEsXLgQa9euxc6dO+tNqp0yZYrt/3v16oX4+HjExMRg586dGDVqlBRRr2vw4MEYPHiw7eMhQ4age/fu+Pjjj/H6669LmKz1LFu2DL169cLAgQPrPe+I7xeRlJTU9gFs/+o40ntmL4rrqWvfvj1cXFxQUFBQ7/mCggKEhIQ0+JqQkJBrXl/336b8mW2hOfdW57333sPChQvx888/Iz4+/prXdurUCe3bt0d6enqLMzdGS+6rjpubGxISEmyZ5fCeteS+ysvLsXbtWjz66KPX/Tr2fr+a62rfZ76+vvDw8GiVfwfOTqntn1LbPoDtX0OU1v7Zs+1TXFGn0WjQr18/bN++3facxWLB9u3b6/1m80eDBw+udz0AbN261XZ9dHQ0QkJC6l1jMBhw4MCBq/6ZbaE59wYA77zzDl5//XVs3rwZ/fv3v+7XycnJQXFxMUJDQ1sl9/U0977+yGw2Izk52ZZZDu9ZS+7rm2++gdFoxIMPPnjdr2Pv96u5rvd91hr/DpydUts/pbZ9ANu/hiit/bNr29ekZRUOYu3atUKr1YoVK1aIkydPiieeeEL4+fnZlnxPmzZNvPTSS7brf/31V+Hq6iree+89cerUKfHqq682uKTfz89PfPfdd+L48eNi8uTJkm1p0pR7W7hwodBoNOLbb7+ttwS8tLRUCCFEaWmpeOGFF8S+fftEZmam2LZtm+jbt6/o3LmzqKqqku19zZ8/X2zZskWcPXtWJCYmiilTpgh3d3dx4sSJevcu9XvW1PuqM2zYMHHfffdd8bxc3q+6LEePHhVHjx4VAMT7778vjh49KrKysoQQQrz00kti2rRptuvrlvW/+OKL4tSpU2Lx4sUNLuu/1t8XXZ9S2z+ltn3NuTe2f9K+Z3Ju+xRZ1AkhxIcffigiIyOFRqMRAwcOFPv377d9bvjw4WL69On1rv/6669Fly5dhEajEXFxceLHH3+s93mLxSJeeeUVERwcLLRarRg1apRITU21x61coSn31rFjRwHgiserr74qhBCioqJCjBkzRgQGBgo3NzfRsWNH8fjjj0vyQ7Qp9/Xcc8/Zrg0ODha33nqrOHLkSL0/Ty7vWVP/LZ4+fVoAED///PMVf5ac3q8dO3Y0+G+r7n6mT58uhg8ffsVr+vTpIzQajejUqVO9/afqXOvvixpHqe2fUts+Idj+1XGE9k/ObZ9KCCGa1rdHRERERHKjuDl1RERERM6IRR0RERGRArCoIyIiIlIAFnVERERECsCijoiIiEgBWNQRERERKQCLOiIiIiIFYFFHREREpAAs6oiIiIgUgEUdOaznn38ed955p9QxiIjsim0fXQ2LOnJYBw8eRP/+/aWOQURkV2z76Gp49is5nOrqanh5eaGmpsb23KBBg7B//34JUxERtS22fXQ9rlIHIGoqV1dX/Prrrxg0aBCSkpIQHBwMd3d3qWMREbUptn10PSzqyOGo1Wrk5eUhICAAvXv3ljoOEZFdsO2j6+GcOnJIR48eZaNGRE6HbR9dC4s6ckhJSUls2IjI6bDto2thUUcOKTk5GX369JE6BhGRXbHto2thUUcOyWKxIDU1FXl5edDr9VLHISKyC7Z9dC0s6sghvfHGG1ixYgU6dOiAN954Q+o4RER2wbaProX71BEREREpAHvqiIiIiBSARR0RERGRArCoIyIiIlIAFnVERERECsCijoiIiEgBWNQRERERKQCLOiIiIiIFYFFHREREpAAs6oiIiIgUgEUdERERkQKwqCMiIiJSABZ1RERERArw/wG4APAEUxDXOQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 2)\n",
    "axs[0].plot(t_eval, solution[\"y squared\"](t_eval))\n",
    "axs[0].set_ylabel(r\"$y^2$\")\n",
    "axs[0].set_xlabel(r\"$t$\")\n",
    "axs[1].plot(solution.t, solution[\"y squared\"].sensitivities[\"a\"])\n",
    "axs[1].set_ylabel(r\"$\\frac{dy^2}{da}$\")\n",
    "axs[1].set_xlabel(r\"$t$\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "495bcf05",
   "metadata": {},
   "source": [
    "## Sensitivities for the DFN model\n",
    "\n",
    "We can do the same for the DFN model included in PyBaMM. We will setup the DFN model using \"Current function\" as an input parameter. This is the parameter we wish to calculate the sensitivities with respect to."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e9119add",
   "metadata": {},
   "outputs": [],
   "source": [
    "# now lets do the same for the DFN model\n",
    "\n",
    "# load model\n",
    "model = pybamm.lithium_ion.DFN()\n",
    "\n",
    "# load parameter values and process model and geometry\n",
    "param = model.default_parameter_values\n",
    "\n",
    "# we want to calculate the sensitivities of the \"Current function\" parameter, so set\n",
    "# this an an input parameter\n",
    "param.update({\"Current function [A]\": \"[input]\"})\n",
    "\n",
    "solver = pybamm.IDAKLUSolver(rtol=1e-3, atol=1e-6)\n",
    "\n",
    "sim = pybamm.Simulation(model, parameter_values=param, solver=solver)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f198fbfe",
   "metadata": {},
   "source": [
    "We can now evaluate the senstivities of, for example, the \"Terminal voltage\" output of the model with respect to the input parameter \"Current function\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "1d794537",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAGwCAYAAAC5ACFFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABsoElEQVR4nO3deViUVf8G8HsGGPZddgFBVERxT8Q1g1yw1PRnoeSShmaapWZpqaW+aWaZaZaZuaVmZVppqeG+IbiACiIKmMgyICIM+zbP7w90ihCcgRmGgftzXXO98mxzH4fX+Xae85wjEgRBABERERE9lljbAYiIiIgaMxZLRERERLVgsURERERUCxZLRERERLVgsURERERUCxZLRERERLVgsURERERUC31tB2gK5HI50tLSYG5uDpFIpO04REREpARBEJCXlwdnZ2eIxTX3H7FYUoO0tDS4urpqOwYRERHVwd27d9GyZcsa97NYUgNzc3MAlX/ZFhYWWk5DREREypDJZHB1dVV8j9eExZIaPLr1ZmFhwWKJiIhIxzxpCA0HeBMRERHVgsUSERERUS1YLBERERHVgsUSERERUS1YLBERERHVgsUSERERUS3qNHWAXC5HQkICMjMzIZfLq+zr37+/WoIRERERNQYqF0vnz5/HuHHjcOfOHQiCUGWfSCRCRUWF2sIRERERaZvKxdJrr72GHj164I8//oCTkxPXQiMiIqImTeVi6datW9izZw+8vLw0kYeIiIioUVF5gLefnx8SEhI0kYWIiIio0VG5Z+mNN97A3LlzIZVK4evrCwMDgyr7O3XqpLZwRERERNomEv47SvsJxOLqnVEikQiCIDTbAd4ymQyWlpbIzc1V60K6DwpKkVNUBidLIxgZ6KntukRERKT897fKPUu3b9+uVzBS3nPrziA1pwj7Xu+Nrm7W2o5DRETULKlcLLm7u2siBz2GpbEBUnOKkFtUpu0oREREzVadJqVMTEzEmjVrEBcXBwDw8fHBm2++idatW6s1XHNnaVw5HiynkMUSERGRtqj8NNzhw4fh4+ODyMhIdOrUCZ06dUJERAQ6dOiAsLAwTWRstlqYGwIAsvJLtJyEiIio+VK5Z2n+/PmYPXs2Pv7442rb3333XTz77LNqC9fc2T8sljJkxVpOQkRE1Hyp3LMUFxeHKVOmVNs+efJkXL9+XS2hqJKLlTEAIDWnSMtJiIiImi+ViyU7OztER0dX2x4dHQ17e3t1ZKKHXG1MAADJ2YVaTkJERNR8qXwbLjQ0FFOnTkVSUhJ69+4NADh79ixWrlyJOXPmqD1gc9bKtrJYupNVqJjHioiIiBqWyj1LixYtwuLFi7Fu3ToMGDAAAwYMwJdffokPP/wQCxcu1ERGAEB2djZCQkJgYWEBKysrTJkyBfn5+bWes3HjRjz99NOwsLCASCRCTk5OtWNatWoFkUhU5fXf8Vja4mZrArEIyCspR2YeB3kTERFpg8rFkkgkwuzZs5GSkoLc3Fzk5uYiJSUFb775pkZ7PkJCQhAbG4uwsDAcOHAAp06dwtSpU2s9p7CwEEOGDMF7771X63FLly5Fenq64vXGG2+oM3qdGerroVULUwBAvDRPy2mIiIiapzrNs/SIubm5unLUKi4uDocOHcKFCxfQo0cPAMC6desQFBSETz/9FM7Ozo8976233gIAnDhxotbrm5ubw9HRUek8JSUlKCn5p6dHJpMpfa6q2jtaIOleAW5IZejf1k5j70NERESPp1TPUrdu3fDgwQMAQNeuXdGtW7caX5oQHh4OKysrRaEEAIGBgRCLxYiIiKj39T/++GPY2tqia9euWLVqFcrLy2s9fsWKFbC0tFS8XF1d652hJt6OlQVpbJrmCjIiIiKqmVI9SyNGjIChoaHizw090FgqlVZ70k5fXx82NjaQSqX1uvasWbPQrVs32NjY4Ny5c1iwYAHS09OxevXqGs9ZsGBBlcHsMplMYwWTb0tLAMC1lFyNXJ+IiIhqp1Sx9MEHHyj+/OGHH6rtzefPn4+VK1fWesyjJVU05d9FT6dOnSCRSDBt2jSsWLFCUSD+l6GhYY371K1zSysAQFJWAbILSmFjKmmQ9yUiIqJKKo9Z8vT0xIULF2Bra1tle05ODrp164akpCSlrzV37lxMmjTpie/n6OiIzMzMKtvLy8uRnZ2t0lgjZfj5+aG8vBx///032rVrp9Zr14W1qQRe9mZIyMzHpTsP8KyPg7YjERERNSsqF0t///03Kioqqm0vKSlBSkqKSteys7ODnd2TBy37+/sjJycHly5dQvfu3QEAx44dg1wuh5+fn0rv+STR0dEQi8WNaoLNHu7WSMjMx8W/s1ksERERNTCli6Xff/9d8efDhw/D0tJS8XNFRQWOHj0KDw8P9aZ7qH379hgyZAhCQ0OxYcMGlJWVYebMmQgODlY8CZeamoqAgABs374dPXv2BFA51kkqlSIhIQEAcO3aNZibm8PNzQ02NjYIDw9HREQEBg4cCHNzc4SHh2P27Nl4+eWXYW1trZG21EUvT1vsvnAXZxKytB2FiIio2VG6WBo5ciSAynmWJk6cWGWfgYEBWrVqhc8++0yt4f5t586dmDlzJgICAiAWizF69GisXbtWsb+srAzx8fEoLPxnaZANGzZgyZIlip/79+8PANiyZQsmTZoEQ0ND7N69Gx9++CFKSkrg4eGB2bNnN7qZyPt4tQBQ+URcVn4JWpg1zHgpIiIiAkSCIAiqnODh4YELFy6gRYsWmsqkc2QyGSwtLZGbmwsLCwuNvEfQF6dxPV2Gz8Z0xujuLTXyHkRERM2Jst/fKs/gffv2bRZKWhD4cKzS4dj6TZVAREREqlG5WJo1a1aV21+PfPnll4oZs0n9BneoLJZO3bqHgpLaJ80kIiIi9VG5WPrll1/Qp0+fatt79+6NPXv2qCUUVefjZAGPFqYoLpPjr+vsXSIiImooKhdL9+/fr/Ik3CMWFhbIyuLTWpoiEokwvHPlk3+/XErVchoiIqLmQ+ViycvLC4cOHaq2/eDBg/D09FRLKHq8/+veEiIRcCYhC7ezCrQdh4iIqFlQeVLKOXPmYObMmbh37x6eeeYZAMDRo0fx2WefYc2aNerOR//iamOCge3scexGJnaev4OFz/loOxIREVGTp3KxNHnyZJSUlOCjjz7CsmXLAACtWrXC119/jQkTJqg9IFX1ci83HLuRiZ8vpWDOoLYwkaj8ERIREZEKVJ5n6d/u3bsHY2NjmJmZqTOTzmmIeZYeqZALGPjpCSRnF2LRcz6Y0lczs6YTERE1dRqbZ+nf7Ozsmn2h1ND0xCJMf7o1AGDjqUSUlFdfp4+IiIjUR+V7OB4eHhCJRDXuT0pKqlcgerJR3Vyw9ugtpOcW46cLdzHev5W2IxERETVZKhdL/514sqysDFFRUTh06BDmzZunrlxUC0N9PUx/ujUW/xaLL44mYFS3ljA15NglIiIiTVD5G/bNN9987Pb169fj4sWL9Q5Eygl+yg3fnbmNO/cLsen0bbwZ2EbbkYiIiJqkeo1Z+rehQ4fil19+Udfl6Akk+mK8PagdgMqxS1n5JVpORERE1DSprVjas2cPbGxs1HU5UsIwXyd0bmmJgtIKrDoUr+04RERETZLKt+G6du1aZYC3IAiQSqW4d+8evvrqK7WGo9qJxSIsfr4DRn99Dj9evIsXn2qJ7u4sWImIiNRJ5WJp5MiRVX4Wi8Wws7PD008/DW9vb3XlIiV1d7fGiz1a4qeLKVj4ayz2z+wDfT21dRgSERE1e0oVS3PmzMGyZctgamqKgQMHwt/fHwYGBprORkqaP7Q9DsdmIC5dho2nk/D6017ajkRERNRkKNUFsW7dOuTn5wMABg4ciAcPHmg0FKnGxlSCRQ/XiVtz5BYSMvO1nIiIiKjpUKpnqVWrVli7di0GDRoEQRAQHh4Oa2vrxx7bv39/tQYk5Yzu5oL9V9Jw8uY9zNtzBT9P8+ftOCIiIjVQam24X3/9Fa+99hoyMzMhEolQ0ykikQgVFc1v+Y2GXBuuNqk5RRjy+SnklZRj3uB2mDGQt+OIiIhqota14UaOHAmpVAqZTAZBEBAfH48HDx5Ue2VnZ6utAaQ6FytjfDC8AwDg87CbiEnN1XIiIiIi3afSfRozMzMcP34cHh4esLS0fOyLtGt0NxcM7uCAcrmAN3dHoai0+fX0ERERqZPKg1oGDBgAfX2uQ9ZYiUQirBjVCfbmhki8V4D//XFd25GIiIh0GkcAN0E2phJ8/lIXiETAzohkHIqRajsSERGRzmKx1ET18WqBqf09AQDv7LmCu9mFWk5ERESkm1gsNWFvD2qHrm5WkBWXY+YPUSgtl2s7EhERkc5RuViaPHky8vLyqm0vKCjA5MmT1RKK1MNAT4x1Y7vC0tgAV+7mYNXhG9qOREREpHNULpa2bduGoqKiatuLioqwfft2tYQi9WlpbYJV/9cJAPDt6ds4Gpeh5URERES6ReliSSaTITc3F4IgIC8vDzKZTPF68OAB/vzzT9jb22syK9XRoA6OeKVPKwDAnJ+uIOUBxy8REREpS+k5AKysrCASiSASidC2bdtq+0UiEZYsWaLWcKQ+84d643JyDq7czcGMXVH4eZo/JPocskZERPQkShdLx48fhyAIeOaZZ/DLL7/AxsZGsU8ikcDd3R3Ozs4aCUn1Z6ivhy/HdsVz687gyt0cLP8zDh8+nO2biIiIaqbU2nD/dufOHbi6ukIsZq/EI41lbThlHI3LwJRtFwEAX4V0Q5Cvk5YTERERaYey398qT8Xt7u6OnJwcREZGIjMzE3J51cfRJ0yYoHpaajAB7R3w2oDW2HAyEe/suQpvR3N42plpOxYREVGjpXLP0v79+xESEoL8/HxYWFhAJBL9czGRqFkupqtLPUsAUF4hx7hNEYi8nQ1vR3Pse70PjCV62o5FRETUoJT9/lb5XtrcuXMxefJk5OfnIycnBw8ePFC8NFkoZWdnIyQkBBYWFrCyssKUKVOQn59f6/FvvPEG2rVrB2NjY7i5uWHWrFnIzc2tclxycjKGDRsGExMT2NvbY968eSgvL9dYOxoDfT0xvhzbFS3MDHFDmof3912DijUzERFRs6FysZSamopZs2bBxMREE3lqFBISgtjYWISFheHAgQM4deoUpk6dWuPxaWlpSEtLw6effoqYmBhs3boVhw4dwpQpUxTHVFRUYNiwYSgtLcW5c+ewbds2bN26FYsXL26IJmmVvYURvhzXFXpiEfZGpWJXZLK2IxERETVKKt+GGzVqFIKDg/Hiiy9qKlM1cXFx8PHxwYULF9CjRw8AwKFDhxAUFISUlBSln8L7+eef8fLLL6OgoAD6+vo4ePAgnnvuOaSlpcHBwQEAsGHDBrz77ru4d+8eJBKJUtfVtdtw/7bhZCI+PngDEj0xfn7NH51drbQdiYiIqEFobID3sGHDMG/ePFy/fh2+vr4wMDCosn/48OGqp32C8PBwWFlZKQolAAgMDIRYLEZERAReeOEFpa7z6C9DX19fcV1fX19FoQQAgwcPxvTp0xEbG4uuXbs+9jolJSUoKSlR/CyTyerSrEZhWn9PXL7zAH9dz8DrOy9j/xt9YWOqXJFIRETUHKhcLIWGhgIAli5dWm2fSCRCRUVF/VP9h1QqrTY7uL6+PmxsbCCVSpW6RlZWFpYtW1bl1p1UKq1SKAFQ/FzbdVesWNFkJuAUiUT49MXOGL7uDP6+X4g3d0dh6ys9oScWPflkIiKiZkDlMUtyubzGl6qF0vz58xWzgtf0unGj/ou/ymQyDBs2DD4+Pvjwww/rfb0FCxYgNzdX8bp79269r6lNFkYG2DC+O4wMxDh9Kwufh93UdiQiIqJGQ+WepX8rLi6GkZFRnc+fO3cuJk2aVOsxnp6ecHR0RGZmZpXt5eXlyM7OhqOjY63n5+XlYciQITA3N8e+ffuq3DZ0dHREZGRkleMzMjIU+2piaGgIQ0PDWt9X13g7WmDl6E54c3c0vjyegM6uVnjWx+HJJxIRETVxKvcsVVRUYNmyZXBxcYGZmRmSkpIAAIsWLcJ3332n0rXs7Ozg7e1d60sikcDf3x85OTm4dOmS4txjx45BLpfDz8+vxuvLZDIMGjQIEokEv//+e7XCzt/fH9euXatSiIWFhcHCwgI+Pj4qtaUpGNHFBZN6twIAzPkxGrezCrQbiIiIqBFQuVj66KOPsHXrVnzyySdVnhbr2LEjNm3apNZwj7Rv3x5DhgxBaGgoIiMjcfbsWcycORPBwcGKJ+FSU1Ph7e2t6Cl6VCgVFBTgu+++g0wmg1QqhVQqVdwuHDRoEHx8fDB+/HhcuXIFhw8fxsKFCzFjxowm13OkrPeC2qOHuzXySsoxfcclFJY27TmniIiInkTlYmn79u3YuHEjQkJCoKf3z6zPnTt3Vsv4oprs3LkT3t7eCAgIQFBQEPr27YuNGzcq9peVlSE+Ph6FhYUAgMuXLyMiIgLXrl2Dl5cXnJycFK9HY4z09PRw4MAB6Onpwd/fHy+//DImTJjw2MHrzYVEX4z1Id0UE1Yu2MsJK4mIqHlTeZ4lY2Nj3LhxA+7u7jA3N8eVK1fg6emJ69evo2fPnrXOqt1U6fI8SzWJSLqPcZsiUCEXsGR4B0x8eHuOiIioqdDYcic+Pj44ffp0te179uypcV4i0j1+nrZYMNQbALDswHVc/Lv5rflHREQE1OFpuMWLF2PixIlITU2FXC7H3r17ER8fj+3bt+PAgQOayEhaMqWvB6Lu5uCPq+mYvvMy/nijL+wt6v70IxERkS5SuWdpxIgR2L9/P44cOQJTU1MsXrwYcXFx2L9/P5599llNZCQtEYlE+GR0J7R1MMO9vBJM33kZpeVybcciIiJqUCoVS+Xl5Vi6dCk8PDwQFhaGzMxMFBYW4syZMxg0aJCmMpIWmRrq45vxPWBupI9Ldx5g2YHr2o5ERETUoFQqlvT19fHJJ5+gvJyPkzcnHi1MsealLgCA78/fwc8XdXvGciIiIlWofBsuICAAJ0+e1EQWasQC2jvgrcA2AID3f43BtZRcLSciIiJqGCoP8B46dCjmz5+Pa9euoXv37jA1Na2yf/jw4WoLR43LrGfaICY1F0fiMvHajkv4fWYf2Jo1z8k7iYio+VB5niWxuObOKJFIpPJiuk1BU5xnqSay4jKM+PIsbmcVoHdrW2yf3BP6eip3UBIREWmdxuZZksvlNb6aY6HU3FgYGeCb8d1hItHDucT7+ORwvLYjERERaZRKxVJZWRn09fURExOjqTykA9o6mOPTMZ0BABtPJWH/lTQtJyIiItIclYolAwMDuLm5sQeJEOTrhNcGtAYAzNtzBdfTZFpOREREpBkq34Z7//338d577yE7m8tfNHfzBrdDvzYtUFwmx9TvL+JBQam2IxEREamdygO8u3btioSEBJSVlcHd3b3a03CXL19Wa0Bd0JwGeP9XTmEpRqw/izv3C9HHyxbbXuGAbyIi0g3Kfn+rPHXAyJEj65OLmhgrEwk2ju+BF746i7MJ9/HxwRtY+JyPtmMRERGpjco9S1Rdc+5ZeuRQTDpe21HZq/j5S53xQteWWk5ERERUO41NHUD0OEM6OuGNZ7wAAPN/ucYZvomIqMlQuVgSi8XQ09Or8UXN1+zAtgjwtkdJuRzTvr+IrPwSbUciIiKqN5XHLO3bt6/Kz2VlZYiKisK2bduwZMkStQUj3SMWi/B5cBeMXH8WSfcK8PrOy9j5qh8MOOCbiIh0mNrGLO3atQs//vgjfvvtN3VcTqdwzFJVCZn5GLn+LPJLyjHR3x1LRnTUdiQiIqJqGnzMUq9evXD06FF1XY50mJe9Gda81AUAsC38Dn66cFe7gYiIiOpBLcVSUVER1q5dCxcXF3VcjpqAQB8HzA5sCwBY+GsMopIfaDkRERFR3ag8Zsna2hoikUjxsyAIyMvLg4mJCXbs2KHWcKTb3njGC7FpufjregZe23EJ+2f2hb2FkbZjERERqUTlMUvbtm2r8rNYLIadnR38/PxgbW2t1nC6gmOWapZfUo4X1p/Frcx8dHOzwg9Te8FQn09NEhGR9in7/c1JKdWAxVLtbmcVYPiXZ5BXXI6xPd2wYpSvtiMRERGpf4D3rVu3MHbsWMhk1VeXz83Nxbhx45CUlFS3tNSkebQwxdqxXSESAT9EJmNnxB1tRyIiIlKa0sXSqlWr4Orq+tjKy9LSEq6urli1apVaw1HTMbCdPeYNbgcA+PD3WFz4O1vLiYiIiJSjdLF08uRJjBkzpsb9L774Io4dO6aWUNQ0TR/QGsN8nVBWIWD6jstIzy3SdiQiIqInUrpYSk5Ohr29fY37W7Rogbt3OZ8O1UwkEmHVmE7wdjRHVn4Jpm6/hKLSCm3HIiIiqpXSxZKlpSUSExNr3J+QkMDBzfREJhJ9fDuhB2xMJbiWmot3frkKPmNARESNmdLFUv/+/bFu3boa969duxb9+vVTSyhq2lxtTPBVSDfoi0XYfyUNX52ouQgnIiLSNqWLpQULFuDgwYP4v//7P0RGRiI3Nxe5ubmIiIjA6NGjcfjwYSxYsECTWakJ6eVpiyUjOgAAVh2Ox1+xUi0nIiIiejyV5lk6cOAAJk+ejPv371fZbmtri02bNmH48OFqD6gLOM9S3S3+LQbbw+/AVKKHX17vDW9H/v0REVHD0NiklEVFRTh06BASEhIgCALatm2LQYMGwcTEpN6hdRWLpborq5BjwneRCE+6D1cbY/w2oy9sTCXajkVERM0AZ/BuQCyW6udBQSlGrD+L5OxC9PK0wfdT/GCgp5Y1nomIiGqk9hm8iTTF2lSCTRN7wMxQH+eTsrFkf6y2IxERESnoTLGUnZ2NkJAQWFhYwMrKClOmTEF+fn6tx7/xxhto164djI2N4ebmhlmzZiE3N7fKcSKRqNpr9+7dmm4O/UdbB3OseakLRCJgx/lkfH+eS6IQEVHjoDPFUkhICGJjYxEWFoYDBw7g1KlTmDp1ao3Hp6WlIS0tDZ9++iliYmKwdetWHDp0CFOmTKl27JYtW5Cenq54jRw5UoMtoZoE+jgolkRZ8nsswhPvP+EMIiIizdOJMUtxcXHw8fHBhQsX0KNHDwDAoUOHEBQUhJSUFDg7Oyt1nZ9//hkvv/wyCgoKoK+vD6CyZ2nfvn0qFUglJSUoKSlR/CyTyeDq6soxS2ogCALe+jEav0WnwdrEAL/N6As32+b78AAREWmOxsYs6enpITMzs9r2+/fvQ09PT9XLKSU8PBxWVlaKQgkAAgMDIRaLERERofR1Hv1lPCqUHpkxYwZatGiBnj17YvPmzU+cUXrFihWwtLRUvFxdXVVrENVIJBJh5ehO6NzSEg8KyxC6/SLyS8q1HYuIiJoxlYulmgqJkpISSCSaeeRbKpVWW5dOX18fNjY2kEqVm8wwKysLy5Ytq3brbunSpfjpp58QFhaG0aNH4/XXX691pnKgcoLOR5Ny5ubmck08NTMy0MM343vA3twQ8Rl5eGt3NOTyRt8BSkRETZT+kw+ptHbtWgCV/+W/adMmmJmZKfZVVFTg1KlT8Pb2VunN58+fj5UrV9Z6TFxcnErXfByZTIZhw4bBx8cHH374YZV9ixYtUvy5a9euKCgowKpVqzBr1qwar2doaAhDQ8N656KaOVoaYeOEHnjxm3AcicvA6rCbePvheCYiIqKGpHSx9PnnnwOo7FnasGFDlVtuEokErVq1woYNG1R687lz52LSpEm1HuPp6QlHR8dqt/7Ky8uRnZ0NR0fHWs/Py8vDkCFDYG5ujn379sHAwKDW4/38/LBs2TKUlJSwINKyLq5WWDnaF7N/vIIvjyegraM5hndWbnwaERGRuihdLN2+fRsAMHDgQOzbtw9WVlb1fnM7OzvY2dk98Th/f3/k5OTg0qVL6N69OwDg2LFjkMvl8PPzq/E8mUyGwYMHw9DQEL///juMjIye+F7R0dGwtrZmodRIvNC1JW5I8/DNySTM+/kKPGxN4dvSUtuxiIioGVFpzFJZWRmSk5ORnp6uqTyP1b59ewwZMgShoaGIjIzE2bNnMXPmTAQHByuehEtNTYW3tzciIyMBVBZKgwYNQkFBAb777jvIZDJIpVJIpVJUVFQAAPbv349NmzYhJiYGCQkJ+Prrr7F8+XK88cYbDdo+qt07g73xjLc9SsrlCN1+ERmyYm1HIiKiZkSlYsnAwADFxdr5otq5cye8vb0REBCAoKAg9O3bFxs3blTsLysrQ3x8PAoLCwEAly9fRkREBK5duwYvLy84OTkpXo8GZBsYGGD9+vXw9/dHly5d8M0332D16tX44IMPtNJGejw9sQhfBHeBl70ZpLJihG6/iKLSCm3HIiKiZkLleZaWL1+OmzdvYtOmTdUewW+uuDZcw0i+X4gR68/gQWEZgnwd8eXYbhCLRdqORUREOkpjC+m+8MILOHr0KMzMzODr6wtTU9Mq+/fu3Vu3xDqMxVLDufB3NkK+jUBphRwzB3rxCTkiIqozZb+/Ve4asrKywujRo+sVjqiunmplgxWjfDH358on5DztTDGqW0ttxyIioiZM5WJpy5YtmshBpLTR3VsiKSsf648nYv4v1+BqY4KnWtloOxYRETVRKs/gvXnzZsU0AkTaMvfZdhja0RGlFXJM+/4Sku8XajsSERE1USoXSytWrICXlxfc3Nwwfvx4bNq0CQkJCZrIRlQjsViE1S92ga+LJbILSjF52wXIisu0HYuIiJoglYulW7duITk5GStWrICJiQk+/fRTtGvXDi1btsTLL7+siYxEj2Us0cOmiT3gaGGEhMx8zNh5GeUVcm3HIiKiJkblp+H+rbCwEKdPn8YPP/yAnTt3QhAElJc3vxXi+TScdsWk5mLMhnAUlVVggr87lo7oqO1IRESkA5T9/la5Z+mvv/7Ce++9h969e8PW1hYLFiyAtbU19uzZg3v37tUrNFFddHSxxJrgLhCJgO3hd7Dt3N/ajkRERE2Iyj1LYrEYdnZ2mDt3LqZOnaqWNeJ0HXuWGodvTiZixcEbEIuAzZOewtPt7LUdiYiIGjGN9SytXr0affr0wSeffIIOHTpg3Lhx2LhxI27evFmvwET1NbW/J17s0RJyAZi5Kwrx0jxtRyIioiagXmOWrl27hpMnT+LYsWM4cOAA7O3tkZKSos58OoE9S41Habkc47+LQMTtbLS0NsavM/qghZmhtmMREVEjpLGeJQAQBAGXL19GWFgYDh8+jOPHj0Mul8POzq7OgYnUQaIvxoaXu6OVrQlSHhRh2veXUFzGRXeJiKjuVC6Wnn/+edja2qJnz57YuXMn2rZti23btiErKwtRUVGayEikEmtTCb6b9BQsjPRx6c4DzP/lKurRgUpERM2cysudeHt7Y9q0aejXrx8sLS01kYmo3lrbmeHrl7tj4uZI/BqdBk87M8wKaKPtWEREpIPqNWaJKnHMUuP1Q2QyFuy9BgBY81IXjOzqouVERETUWGh0zBKRrhjb0w1T+3sCAObtuYLwxPtaTkRERLqGxRI1efOHeCPI1xFlFQKmfX8RCZmcUoCIiJTHYomavEeL7nZzs4KsuByTtlzAvbwSbcciIiIdwWKJmgUjAz1smviUYkqBKdsuoLC0+a1jSEREqqtTsZSYmIiFCxdi7NixyMzMBAAcPHgQsbGxag1HpE42phJsfaUnrE0McDUlF7N+iEKFnM83EBFR7VQulk6ePAlfX19ERERg7969yM/PBwBcuXIFH3zwgdoDEqlTqxam2DSxByT6YhyJy8TS/bGcg4mIiGqlcrE0f/58/O9//0NYWBgkEoli+zPPPIPz58+rNRyRJnR3t8Gal7oAALaF38F3Z25rNxARETVqKhdL165dwwsvvFBtu729PbKystQSikjTgnyd8F6QNwDgoz/jcPBaupYTERFRY6VysWRlZYX09OpfLFFRUXBx4YR/pDtC+3lifC93CALw1o/RuHTngbYjERFRI6RysRQcHIx3330XUqkUIpEIcrkcZ8+exdtvv40JEyZoIiORRohEInzwvA8CvO1RUi5H6PaL+DurQNuxiIiokVG5WFq+fDm8vb3h6uqK/Px8+Pj4oH///ujduzcWLlyoiYxEGqOvJ8a6cV3h62KJ7IJSvLL1Ah4UlGo7FhERNSJ1XhsuOTkZMTExyM/PR9euXdGmTfNdpJRrw+m+TFkxXvjqHFJzitDD3Ro7XvWDkYGetmMREZEGKfv9zYV01YDFUtNwMyMPo78+h7zicgzr5IR1wV0hFou0HYuIiDRE2e9vfVUvPGfOnMduF4lEMDIygpeXF0aMGAEbGxtVL02kVW0dzPHN+O6YuDkSf1xNR0trYywY2l7bsYiISMtU7lkaOHAgLl++jIqKCrRr1w4AcPPmTejp6cHb2xvx8fEQiUQ4c+YMfHx8NBK6sWHPUtOy93IK5vx0BQCwbGRHjO/lruVERESkCcp+f6s8wHvEiBEIDAxEWloaLl26hEuXLiElJQXPPvssxo4di9TUVPTv3x+zZ8+uVwOItGVUt5aY82xbAMDi32LwJ+dgIiJq1lTuWXJxcUFYWFi1XqPY2FgMGjQIqampuHz5MgYNGtRsJqlkz1LTIwgC3v81BrsikiHRE2PzpKfQt00LbcciIiI10ljPUm5urmLx3H+7d+8eZDIZgMqJK0tL+fg16S6RSIRlIzoiyNcRpRVyTP3+Iq7czdF2LCIi0oI63YabPHky9u3bh5SUFKSkpGDfvn2YMmUKRo4cCQCIjIxE27Zt1Z2VqEHpiUX4/KUu6OvVAoWlFZi0JRIJmfnajkVERA1M5WLpm2++QUBAAIKDg+Hu7g53d3cEBwcjICAAGzZsAAB4e3tj06ZNag2anZ2NkJAQWFhYwMrKClOmTEF+fu1fXNOmTUPr1q1hbGwMOzs7jBgxAjdu3KhyTHJyMoYNGwYTExPY29tj3rx5KC8vV2t20l2G+nrYML47Ore0xIPCMkz4LgJpOUXajkVERA2ozvMs5efnIykpCQDg6ekJMzMztQb7r6FDhyI9PR3ffPMNysrK8Morr+Cpp57Crl27ajxn48aN8Pb2hpubG7Kzs/Hhhx8iOjoat2/fhp6eHioqKtClSxc4Ojpi1apVSE9Px4QJExAaGorly5crnY1jlpq+7IJSjNlwDon3CtDazhQ/v9YbNqYSbcciIqJ6aFKTUsbFxcHHxwcXLlxAjx49AACHDh1CUFAQUlJS4OzsrNR1rl69is6dOyMhIQGtW7fGwYMH8dxzzyEtLQ0ODg4AgA0bNuDdd9/FvXv3IJEo92XIYql5SM0pwv99fQ7pucXo3NISO0N7wcxQ5anKiIiokdDYAG8AuHjxIt555x0EBwdj1KhRVV6aEB4eDisrK0WhBACBgYEQi8WIiIhQ6hoFBQXYsmULPDw84Orqqriur6+volACgMGDB0MmkyE2NrbGa5WUlEAmk1V5UdPnYmWM76f0hLWJAa6k5OK17y+hpLxC27GIiEjDVC6Wdu/ejd69eyMuLg779u1DWVkZYmNjcezYMVhaWmoiI6RSKezt7ats09fXh42NDaRSaa3nfvXVVzAzM4OZmRkOHjyIsLAwRY+RVCqtUigBUPxc23VXrFgBS0tLxetR8UVNn5e9Oba80hMmEj2cScjCnB+voELe6DtniYioHlQulpYvX47PP/8c+/fvh0QiwRdffIEbN27gxRdfhJubm0rXmj9/PkQiUa2v/w7IVlVISAiioqJw8uRJtG3bFi+++CKKi4vrdc0FCxYgNzdX8bp79269rke6pYurFTaO7wEDPRH+uJaORb/FQAfuZhMRUR2pPOAiMTERw4YNAwBIJBIUFBRAJBJh9uzZeOaZZ7BkyRKlrzV37lxMmjSp1mM8PT3h6OhYbW6n8vJyZGdnw9HRsdbzH/X+tGnTBr169YK1tTX27duHsWPHwtHREZGRkVWOz8jIAIBar2toaAhDQ8Na35eatr5tWmDNS10x84fL2BWRDFtTCeYOaqftWEREpAEqF0vW1tbIy8sDUDmbd0xMDHx9fZGTk4PCwkKVrmVnZwc7O7snHufv74+cnBxcunQJ3bt3BwAcO3YMcrkcfn5+Sr+fIAgQBAElJSWK63700UfIzMxU3OYLCwuDhYVFs1nXjupuWCcn5BR1xPv7YrDuWAJsTCV4pY+HtmMREZGaqXwbrn///ggLCwMAjBkzBm+++SZCQ0MxduxYBAQEqD0gALRv3x5DhgxBaGgoIiMjcfbsWcycORPBwcGKJ+FSU1Ph7e2t6ClKSkrCihUrcOnSJSQnJ+PcuXMYM2YMjI2NERQUBAAYNGgQfHx8MH78eFy5cgWHDx/GwoULMWPGDPYckVJC/Nwx9+E6ckv2X8evUalaTkREROqmcs/Sl19+qRjz8/7778PAwADnzp3D6NGjsXDhQrUHfGTnzp2YOXMmAgICIBaLMXr0aKxdu1axv6ysDPHx8YreLSMjI5w+fRpr1qzBgwcP4ODggP79++PcuXOKXiQ9PT0cOHAA06dPh7+/P0xNTTFx4kQsXbpUY+2gpmfmM17ILizFlrN/4+2fr8DS2AADve2ffCIREekEnZhnqbHjPEsklwuY+/MV7ItKhZGBGDum+KFHKxttxyIiolpobJ4lPT29xy6ke//+fejp6al6OaImQSwW4ZP/64SB7exQXCbH5K0XcEPK+beIiJoClYulmjqiSkpKlJ7xmqgpMtAT46uQ7ujhbg1ZcTkmfBeJu9mqPfRARESNj9Jjlh6NDxKJRNi0aVOVteAqKipw6tQpeHt7qz8hkQ4xlujhu4lP4aWN4bghzcPL30Vgz2u9YWfOBwaIiHSV0mOWPDwqH4m+c+cOWrZsWeWWm0QiQatWrbB06VKVHuVvKjhmif4rQ1aM/9twDnezi+DjZIHd03rBwshA27GIiOhfNLaQ7sCBA7F3715YW1vXO2RTwWKJHufvrAL834ZwZOWXoKeHDbZP7gkjA47rIyJqLDQ2wPv48eMslIiU0KqFKbZNfgrmhvqIvJ2NmbuiUF4h13YsIiJSkVI9S3PmzFH6gqtXr65XIF3EniWqTUTSfUzYHImScjn+r3tLrPq/ThCJRNqORUTU7Cn7/a3UAO+oqCil3pRfAETV+Xna4stx3fDajkvYcykFNqYSvBfUXtuxiIhISZyUUg3Ys0TK+PniXczbcxUAMH+oN14b0FrLiYiImjeNjVn6t5SUFKSkpNTnEkTNxpgernj/YY/Sxwdv4McLyVpOREREylC5WJLL5Vi6dCksLS3h7u4Od3d3WFlZYdmyZZDLOXiVqDah/T0VPUoL9l7DoRiplhMREdGTqFwsvf/++/jyyy/x8ccfIyoqClFRUVi+fDnWrVuHRYsWaSIjUZPy7pB2eKmHK+QCMGt3FM4lZmk7EhER1ULlMUvOzs7YsGEDhg8fXmX7b7/9htdffx2pqalqDagLOGaJVFVeIceMXZdxODYDZob62D21Fzq6WGo7FhFRs6KxMUvZ2dmPXdbE29sb2dnZql6OqFnS1xPji+Cu8Pe0RX5JOSZujkTSvXxtxyIiosdQuVjq3Lkzvvzyy2rbv/zyS3Tu3FktoYiaAyMDPWyc0B0dXSxwv6AU47+LRHpukbZjERHRf6h8G+7kyZMYNmwY3Nzc4O/vDwAIDw/H3bt38eeff6Jfv34aCdqY8TYc1UdWfgnGbAjH7awCuNua4IfQXnC2MtZ2LCKiJk9jt+EGDBiA+Ph4vPDCC8jJyUFOTg5GjRqF+Pj4ZlkoEdVXCzND7HjVD642xrhzvxDBG88jLYc9TEREjQUnpVQD9iyROqTmFGHsxvNIzi6Em40JfpjaCy7sYSIi0hiN9Sx5eXnhww8/xK1bt+oVkIiqcrEyxu6pveBua4Lk7EIEbwxHKnuYiIi0TuViacaMGfjjjz/Qrl07PPXUU/jiiy8glXJiPSJ1cP5XwXQ3uwjBG8OR8qBQ27GIiJo1lYul2bNn48KFC7hx4waCgoKwfv16uLq6YtCgQdi+fbsmMhI1K06WlQVTK0XBdB53s1kwERFpi1rGLJ0/fx7Tp0/H1atXUVFRoY5cOoVjlkgTpLnFCN4Yjr/vFypu0bnamGg7FhFRk9EgC+lGRkbirbfewgsvvICbN29izJgx9bkcEf2Lo6URdk/1h0cLU6TmsIeJiEhbVC6Wbt68iQ8++ABt27ZFnz59EBcXh5UrVyIjIwO7d+/WREaiZsvR0gg/hPZiwUREpEUq34YTi8V46qmnMG7cOAQHB8PBwUFT2XQGb8ORpmXIijF243kkZRXAxcoYP4T2gpstb8kREdWHst/fKhdLt27dQps2beodsClhsUQN4d8Fk/PDW3QsmIiI6k5jY5ZYKBFph4OFEXZP7QVPO1OkPRz8fed+gbZjERE1efUa4E1EDcvewgi7Q3uhtaJgOs+CiYhIw1gsEekYewsj/DC1smBKzy3GS9+cx99ZLJiIiDSFxRKRDrI3rxyz5GVvBqmssofpNgsmIiKNYLFEpKPszA3xQ2gvtFEUTOEsmIiINECpp+HmzJmj9AVXr15dr0C6iE/DkTbdyyvBuG/P41ZmPhwsKgsoTzszbcciImr0lP3+1lfmYlFRUUq9qUgkUi4dEamNnbkhfpjaC+O+PY+bGfkI3nj+4VNzLJiIiNRBLWvDNXfsWaLGICu/BCHfRiA+Iw/2Dwuo1iyYiIhq1CBrwxFR49HCzBC7Qv3QzsEcmXklGLvxPBLv5Ws7FhGRzqtTsXTx4kW88847CA4OxqhRo6q8NCU7OxshISGwsLCAlZUVpkyZgvz82r8Ipk2bhtatW8PY2Bh2dnYYMWIEbty4UeUYkUhU7cU17khX2T4smLwdKwum4I3nkZDJgomIqD5ULpZ2796N3r17Iy4uDvv27UNZWRliY2Nx7NgxWFpaaiIjACAkJASxsbEICwvDgQMHcOrUKUydOrXWc7p3744tW7YgLi4Ohw8fhiAIGDRoECoqKqoct2XLFqSnpyteI0eO1Fg7iDTN1swQO1+tLJjuKQqmPG3HIiLSWSqPWerUqROmTZuGGTNmwNzcHFeuXIGHhwemTZsGJycnLFmyRO0h4+Li4OPjgwsXLqBHjx4AgEOHDiEoKAgpKSlwdnZW6jpXr15F586dkZCQgNatWwOo7Fnat2+fSgVSSUkJSkpKFD/LZDK4urpyzBI1KtkFpRj37XnckOahhZkhfgj1QxsHc23HIiJqNDQ2ZikxMRHDhg0DAEgkEhQUFEAkEmH27NnYuHFj3RPXIjw8HFZWVopCCQACAwMhFosRERGh1DUKCgqwZcsWeHh4wNXVtcq+GTNmoEWLFujZsyc2b96MJ9WPK1asgKWlpeL13+sRNQY2phLsCu2F9k4WyMovwdhvz+NWBnuYiIhUpXKxZG1tjby8yn9wXVxcEBMTAwDIyclBYWGhetM9JJVKYW9vX2Wbvr4+bGxsIJVKaz33q6++gpmZGczMzHDw4EGEhYVBIpEo9i9duhQ//fQTwsLCMHr0aLz++utYt25drddcsGABcnNzFa+7d+/WvXFEGmRjKsGuV/3g42SBrPxSjP32PG6yYCIiUonKxVL//v0RFhYGABgzZgzefPNNhIaGYuzYsQgICFDpWvPnz3/sAOt/v/47IFtVISEhiIqKwsmTJ9G2bVu8+OKLKC4uVuxftGgR+vTpg65du+Ldd9/FO++8g1WrVtV6TUNDQ1hYWFR5ETVW1qYS7HzVDx2cHxZMG1kwERGpQuUxS9nZ2SguLoazszPkcjk++eQTnDt3Dm3atMHChQthbW2t9LXu3buH+/fv13qMp6cnduzYgblz5+LBgweK7eXl5TAyMsLPP/+MF154Qan3Ky0thbW1NTZt2oSxY8c+9pg//vgDzz33HIqLi2FoaKjUdTnPEumCnMJShGyKQGyaDLYPb9G1c+QYJiJqvtQ6g/e/2djYKP4sFosxf/78uiUEYGdnBzs7uyce5+/vj5ycHFy6dAndu3cHABw7dgxyuRx+fn5Kv58gCBAEocrg7P+Kjo6GtbW10oUSka6wMqnsYXr5uwjEpMow9tvz2DzpKXRxtdJ2NCKiRk3lYgkA5HI5EhISkJmZCblcXmVf//791RLs39q3b48hQ4YgNDQUGzZsQFlZGWbOnIng4GDFk3CpqakICAjA9u3b0bNnTyQlJeHHH3/EoEGDYGdnh5SUFHz88ccwNjZGUFAQAGD//v3IyMhAr169YGRkhLCwMCxfvhxvv/222ttA1BhYmUiwY4ofJmyOxNWUXIzdeB5fhXTDQG/7J59MRNRMqVwsnT9/HuPGjcOdO3eqPTUmEomqzWGkLjt37sTMmTMREBAAsViM0aNHY+3atYr9ZWVliI+PVwwyNzIywunTp7FmzRo8ePAADg4O6N+/P86dO6cYLG5gYID169dj9uzZEAQBXl5eWL16NUJDQzXSBqLGwMqk8hbc9B2XcPpWFl7dfhEfj/LFmB58qpOI6HFUHrPUpUsXtG3bFkuWLIGTk1O1xXM1OTFlY8UxS6SLSsvlmP/LVeyNSgUAvD2oLWYM9OKC2ETUbGhszNKtW7ewZ88eeHl51SsgEWmXRF+Mz17sDHsLI2w4mYhP/7qJDFkJPhzeAXpiFkxERI+oPHWAn58fEhISNJGFiBqYSCTC/KHe+OB5H4hEwPfn72DGzssoLtPM7XQiIl2kcs/SG2+8gblz50IqlcLX1xcGBgZV9nfq1Elt4YioYbzSxwN25oaY8+MVHIqVYsJ3kfh2Qg9Ymhg8+WQioiZO5TFLYnH1ziiRSARBEDQ6wLsx45glairCE+9j6vaLyCspR1sHM2x9pSecrYy1HYuISCOU/f5WuVi6c+dOrfvd3d1VuVyTwGKJmpK4dBkmbo5EZl4JnCyNsG1yT7TlArxE1ARprFii6lgsUVOT8qAQEzdHIvFeASyM9LFp4lPo6WHz5BOJiHSIWoul33//HUOHDoWBgQF+//33Wo8dPny46ml1HIslaooeFJRiyrYLuJycA4m+GGuDu2BIRydtxyIiUhu1FktisRhSqRT29vaPHbOkuBjHLLFYoialqLQCb/wQhSNxGRCJgKXDO2C8fyttxyIiUgtlv7+VmjpALpcrZr2Wy+U1vppjoUTUlBlL9LDh5W4Y29MNggAs+i0Wqw7fqDZ7PxFRU6byPEtE1Lzo64mx/IWOmB3YFgCw/ngi3tlzFWUV8iecSUTUNNRpId0LFy7g+PHjj11Id/Xq1WoJRkSNh0gkwpuBbWBvYYj3913Dz5dScC+/BF+FdIOJpE7/jBAR6QyV/5Vbvnw5Fi5ciHbt2sHBwaHKOlJcU4qoaRvb0w12ZoaY+cNlnIi/h7HfRmDzxB6wNTPUdjQiIo1ReeoABwcHrFy5EpMmTdJQJN3DAd7U3Fy68wBTtl1ATmEZPFqYYtsrPeFma6LtWEREKlHrAO8qJ4jF6NOnT73CEZFu6+5ujT2v9YaLlTFuZxVg1NfnEJOaq+1YREQaoXKxNHv2bKxfv14TWYhIh3jZm2Hv673h7WiOrPwSvPRNOE7fuqftWEREaqfybTi5XI5hw4bh5s2b8PHxqbaQ7t69e9UaUBfwNhw1Z7LiMkzbfgnhSfehLxbh0zGdMbKri7ZjERE9kcZuw82aNQvHjx9H27ZtYWtrC0tLyyovImpeLIwMsHXyUxjWyQnlcgFv/RiNb08laTsWEZHaqNyzZG5ujt27d2PYsGGayqRz2LNEBMjlAv73Rxw2n70NAJjS1wPvB7WHWMynZImocdJYz5KNjQ1at25dr3BE1PSIxSIseq493gvyBgB8d+Y23vwxGiXlnNmfiHSbysXShx9+iA8++ACFhYWayENEOkwkEmFq/9b4/KXO0BeLsP9KGl7ZcgF5xWXajkZEVGcq34br2rUrEhMTIQgCWrVqVW2A9+XLl9UaUBfwNhxRdadu3sP0HZdQUFqB9k4W2PrKU3CwMNJ2LCIiBWW/v1WewXvkyJH1yUVEzUT/tnb4cZo/Jm2JRFy6DMO/PIOvQrqju7u1tqMREalEpWKpvLwcIpEIkydPRsuWLTWViYiaiI4ulvhlem9M2XYRCZn5CN4YjiXDO2Kcn5u2oxERKU2lMUv6+vpYtWoVysvLNZWHiJoYd1tT/DqjD4Z0cERZhYD39l3D/F+ucuA3EekMlQd4P/PMMzh58qQmshBRE2VmqI+vX+6GeYPbQSQCdl+4i5e+OY/03CJtRyMieiKVxywNHToU8+fPx7Vr19C9e3eYmppW2T98+HC1hSOipkMkEmHGQC90cLbArB+iEH03B8+vO4P147rBz9NW2/GIiGqk8tNwYnHNnVEikQgVFc2va51PwxGp5s79Akz7/hJuSPOgLxZh4bD2mNi7FUQiTmBJRA1HY5NSyuXyGl/NsVAiItW525pi7+u98XxnZ5TLBXy4/zrm/nwFxWX8N4SIGh+Vi6V/Ky4uVlcOImpmTCT6WBvcpXJJFBGw93Iq/m/DOaQ84IS3RNS4qFwsVVRUYNmyZXBxcYGZmRmSkioXzFy0aBG+++47tQckoqZLJBIhtL8ndkzxg42pBDGpMjy/7gzOJWRpOxoRkYLKxdJHH32ErVu34pNPPoFEIlFs79ixIzZt2qTWcETUPPT2aoHfZ/ZBRxcLPCgsw8vfReDbU0lQcUglEZFGqFwsbd++HRs3bkRISAj09PQU2zt37owbN26oNRwRNR8trU2w57XeGNXNBXIB+OjPOMzaHY3CUs7rRkTapXKxlJqaCi8vr2rb5XI5ysq4WCYR1Z2RgR4+G9MZS4Z3UCzEO+qrc7hzv0Db0YioGVO5WPLx8cHp06erbd+zZw+6du2qllBE1HyJRCJM7N0Ku0J7oYWZBDekeXh+3RmciM/UdjQiaqaULpYmT56MvLw8LF68GDNnzsTKlSshl8uxd+9ehIaG4qOPPsLixYs1FjQ7OxshISGwsLCAlZUVpkyZgvz8fKXOFQQBQ4cOhUgkwq+//lplX3JyMoYNGwYTExPY29tj3rx5XM6FqBHo6WGDA2/0QxdXK8iKy/HK1gtYfzyB45iIqMEpXSxt27YNRUVFGDFiBPbv348jR47A1NQUixcvRlxcHPbv349nn31WY0FDQkIQGxuLsLAwHDhwAKdOncLUqVOVOnfNmjWPneyuoqICw4YNQ2lpKc6dO4dt27Zh69atGi36iEh5jpZG+HFaL4zt6QpBAFYdjsf0HZeRX8L/oCGihqP0DN5isRhSqRT29vaazlRNXFwcfHx8cOHCBfTo0QMAcOjQIQQFBSElJQXOzs41nhsdHY3nnnsOFy9ehJOTE/bt24eRI0cCAA4ePIjnnnsOaWlpcHBwAABs2LAB7777Lu7du1flab/acAZvIs3bFZGMD36PQVmFAC97M3wzvjta25lpOxYR6TCNzOCdl5cHmUxW60sTwsPDYWVlpSiUACAwMBBisRgRERE1nldYWIhx48Zh/fr1cHR0fOx1fX19FYUSAAwePBgymQyxsbE1XrekpKRB2k1E/xjn54bdU/3hYGGIhMx8jPzyLI5cz9B2LCJqBlQqltq2bQtra+vHvqysrGBtba2RkI/r0dLX14eNjQ2kUmmN582ePRu9e/fGiBEjarzuvwslAIqfa7vuihUrYGlpqXi5uroq2xQiqofu7tbY/0ZfPNXKGnkl5Xh1+0V8HnYTcjnHMRGR5uircvCePXtgY2OjtjefP38+Vq5cWesxcXFxdbr277//jmPHjiEqKqpO59dmwYIFmDNnjuJnmUzGgomogdibG2Hnq73wvz+uY3v4HXxx9BZiUnPxeXAXWBgZaDseETVBKhVLffr0UeuYpblz52LSpEm1HuPp6QlHR0dkZlZ9bLi8vBzZ2dmPvb0GAMeOHUNiYiKsrKyqbB89ejT69euHEydOwNHREZGRkVX2Z2RUduvXdF0AMDQ0hKGhYa25iUhzJPpiLB3REb4ulnj/1xgcvZGJEV+excbx3dHGwVzb8YioiVGpWFI3Ozs72NnZPfE4f39/5OTk4NKlS+jevTuAymJILpfDz8/vsefMnz8fr776apVtvr6++Pzzz/H8888rrvvRRx8hMzNTUQSGhYXBwsICPj4+9WkaETWAMT1c0c7RHK99fwm3swowcv1ZfDqmM4b6Omk7GhE1IUqPWXJ3d6+yvElDat++PYYMGYLQ0FBERkbi7NmzmDlzJoKDgxVPwqWmpsLb21vRU+To6IiOHTtWeQGAm5sbPDw8AACDBg2Cj48Pxo8fjytXruDw4cNYuHAhZsyYwZ4jIh3RqaUV9r/RF708bVBQWoHpOy9j5aEbKK+QazsaETURShdLt2/fhq2trSaz1Grnzp3w9vZGQEAAgoKC0LdvX2zcuFGxv6ysDPHx8SgsLFT6mnp6ejhw4AD09PTg7++Pl19+GRMmTMDSpUs10QQi0hBbM0PsmOKHV/tW/ofQ1ycSMWL9WcSk5mo5GRE1BUrPs0Q14zxLRI3H/itpWPhrDHKLyqAnFuHVvh54K7AtjCXa6RknosZLI/MsERE1ds93dsaROQPwXCcnVMgFfHMqCYPXnMLZhCxtRyMiHcViiYiaHDtzQ3w5rhs2TegBJ0sjJGcXImRTBOb9fAU5haXajkdEOkapYsnGxgZZWZX/VfZoQV0iosYu0McBf83ujwn+7hCJgJ8vpSBw9Sn8cTWdC/ISkdKUKpZKS0sVS3ps27YNxcXFGg1FRKQu5kYGWDqiI36e5g8vezNk5Zdgxq7LCN1+Cem5RdqOR0Q6QKkB3s8++ywyMjLQvXt3bNu2DS+99BKMjY0fe+zmzZvVHrKx4wBvIt1QUl6Br44n4qsTCSirEGBmqI93h3ojpKcbxGKRtuMRUQNT6wDvHTt2ICgoCPn5+RCJRMjNzcWDBw8e+yIiaqwM9fUw+9m2+GNWP3R1s0J+STkW/RqDF78JR0JmvrbjEVEjpfLUAR4eHrh48aJW51xqbNizRKR7KuQCdpy/g08O3UBBaQUkemLMfMYLrw1oDYk+n30hag6U/f7mPEtqwGKJSHel5hRh4b5rOB5/DwDQ1sEMH4/uhG5u1lpORkSaptF5lk6ePInnn38eXl5e8PLywvDhw3H69Ok6hyUi0hYXK2NsnvQUvgjuAhtTCW5m5GP01+ewZH8sCkrKtR2PiBoBlYulHTt2IDAwECYmJpg1axZmzZoFY2NjBAQEYNeuXZrISESkUSKRCCO6uODInAEY1c0FggBsOfs3Bn1+CsfjM7Udj4i0TOXbcO3bt8fUqVMxe/bsKttXr16Nb7/9FnFxcWoNqAt4G46oaTl18x7e23cNKQ8qpxYY2cUZi57zga0ZF9gmako0NmbJ0NAQsbGx8PLyqrI9ISEBHTt2bJZzMLFYImp6CkvLsfqvm9h89jbkAmBtYoDFz/tgZBcXiEScZoCoKdDYmCVXV1ccPXq02vYjR47A1dVV1csRETVKJhJ9LHzOB/te7wNvR3M8KCzD7B+vYOKWC7ibXajteETUgPRVPWHu3LmYNWsWoqOj0bt3bwDA2bNnsXXrVnzxxRdqD0hEpE2dXa2w/42+2HgqCV8cvYVTN+9h0Oen8PbgdpjUuxX0OJklUZNXp6kD9u3bh88++0wxPql9+/aYN28eRowYofaAuoC34Yiah8R7+Viw9xoib2cDqCykVo72hbcj/39PpIs4z1IDYrFE1HzI5QJ2X7iLFX/GIa+kHPpiEaY/3RozBnrByEBP2/GISAUanWeJiKi5EotFGOfnhiNzB2CQjwPK5QLWHUtA0NrTih4nImpaWCwREdWBg4URNk7ogQ0vd4OduSGS7hXgxW/C8f6+a5AVl2k7HhGpEYslIqJ6GNLRCUfmDMDYnpVPA++MSMag1afwV6xUy8mISF1YLBER1ZOlsQFWjOqEXaF+aGVrAqmsGFO/v4TXd15CZl7zm3uOqKmpd7FUUVGB6OhoPHjwQB15iIh0Vu/WLXDorf6Y/nRr6IlF+POaFIGfncRPF+6Cz9IQ6S6Vi6W33noL3333HYDKQmnAgAHo1q0bXF1dceLECXXnIyLSKUYGenh3iDd+n9kHvi6WkBWX451friJkUwT+zirQdjwiqgOVi6U9e/agc+fOAID9+/fj9u3buHHjBmbPno33339f7QGJiHRRB2dL7Hu9N94Pag8jAzHOJd7H4DWnsOFkIsor5NqOR0QqULlYysrKgqOjIwDgzz//xJgxY9C2bVtMnjwZ165dU3tAIiJdpa8nRmh/T/z11gD09WqBknI5Pj54AyPWn0VMaq624xGRklQulhwcHHD9+nVUVFTg0KFDePbZZwEAhYWF0NPjhGxERP/lZmuC76f0xKdjOsPS2ACxaTKMWH8WK/6MQ1FphbbjEdETqFwsvfLKK3jxxRfRsWNHiEQiBAYGAgAiIiLg7e2t9oBERE2BSCTC/3VviSNzBuD5zs6okAv45lQShnxxCucSsrQdj4hqUaflTvbs2YO7d+9izJgxaNmyJQBg27ZtsLKyapbrw3G5EyJS1dG4DCz8NQbpuZVTC7zYoyXeD/KBpYmBlpMRNR8NsjZccXExjIyM6np6k8FiiYjqIq+4DJ8ejsf283cgCEALM0MsGd4BQb6OEIlE2o5H1ORpbG24iooKLFu2DC4uLjAzM0NSUhIAYNGiRYopBYiI6MnMjQywZERH7HnNH172ZsjKL8GMXZcRuv0S0nOLtB2PiB5SuVj66KOPsHXrVnzyySeQSCSK7R07dsSmTZvUGo6IqDno7m6DP2b1xZsBbWCgJ8KRuAw8u/oUvj9/B3I5J7Mk0jaVi6Xt27dj48aNCAkJqfL0W+fOnXHjxg21hiMiai4M9fUw+9m2+GNWP3R1s0J+STkW/RqDF78JR0JmvrbjETVrKhdLqamp8PLyqrZdLpejrIwrbRMR1UdbB3Psea03lgzvAFOJHi7eeYCgL05j3dFbKC3nZJZE2qByseTj44PTp09X275nzx507dpVLaGIiJozPbEIE3u3wl9zBmBgOzuUVsjxWdhNPL/uDKKSuQ4nUUPTV/WExYsXY+LEiUhNTYVcLsfevXsRHx+P7du348CBA5rISETULLlYGWPzpKfw+5U0LN1/HfEZeRj19TkM7+yM0H6e6Ohiqe2IRM2Cyj1LI0aMwP79+3HkyBGYmppi8eLFiIuLw/79+xWzeWtCdnY2QkJCYGFhASsrK0yZMgX5+crdxxcEAUOHDoVIJMKvv/5aZZ9IJKr22r17twZaQESkOpFIhBFdXHBkzgCM7tYSggD8Fp2G59adQcim8zh58x7qMQMMESlB5Z4lAOjXrx/CwsLUnaVWISEhSE9PR1hYGMrKyvDKK69g6tSp2LVr1xPPXbNmTa1zlmzZsgVDhgxR/GxlZaWOyEREamNtKsFnL3bGK31a4dvTSThwNR1nE+7jbMJ9eDuaI7SfJ57v7AyJvsr/DUxET1CvSSkbSlxcHHx8fHDhwgX06NEDAHDo0CEEBQUhJSUFzs7ONZ4bHR2N5557DhcvXoSTkxP27duHkSNHKvaLRKJq21TFSSmJqKGl5hRh85nb2B2ZjIKH68s5WhjhlT6tMNbPDRZGnAmc6EnUOimljY0NsrIq1y6ytraGjY1NjS9NCA8Ph5WVlaJQAoDAwECIxWJERETUeF5hYSHGjRuH9evXw9HRscbjZsyYgRYtWqBnz57YvHnzE7u0S0pKIJPJqryIiBqSi5UxFj3ng3MLAvDuEG/YmxtCKivGioM30HvFMXz0x3Wk5XBiSyJ1UOo23Oeffw5zc3PFnxt6Gn6pVAp7e/sq2/T19WFjYwOpVFrjebNnz0bv3r1rXa9u6dKleOaZZ2BiYoK//voLr7/+OvLz8zFr1qwaz1mxYgWWLFmiekOIiNTM0tgA059ujcl9W+H36DR8ezoJNzPy8e3p29hy9m88/3AwuI8ze72J6kqrt+Hmz5+PlStX1npMXFwc9u7di23btiE+Pr7KPnt7eyxZsgTTp0+vdt7vv/+OuXPnIioqCmZmZgCUu+W2ePFibNmyBXfv3q3xmJKSEpSUlCh+lslkcHV15W04ItI6QRBw4uY9bDyZhPCk+4rt/dq0wNT+nujr1YLrzhE9pOxtOJUHeOvp6SE9Pb1aT8/9+/dhb2+PiooKpa81d+5cTJo0qdZjPD094ejoiMzMzCrby8vLkZ2dXePttWPHjiExMbHaYO3Ro0ejX79+OHHixGPP8/Pzw7Jly1BSUgJDQ8PHHmNoaFjjPiIibRKJRBjYzh4D29njWkouNp5Owp/X0nH6VhZO38pCeycLTO3vgec6OcNAj4PBiZShcs+SWCx+7G2xtLQ0tG7dGkVF6r9H/miA98WLF9G9e3cAwF9//YUhQ4bUOMBbKpUqxlk94uvriy+++ALPP/88PDw8HvteH330ET777DNkZ2crnY8DvImoMbubXYjNZ2/jxwt3UfhwMLiTpREm9/FAcE9XmHMwODVTau9ZWrt2LYDK/2rZtGmT4tYWAFRUVODUqVPw9vauR+SatW/fHkOGDEFoaCg2bNiAsrIyzJw5E8HBwYpCKTU1FQEBAdi+fTt69uwJR0fHx/Y6ubm5KQql/fv3IyMjA7169YKRkRHCwsKwfPlyvP322xppBxGRNrjamOCD5zvgzYA22BmRjK3n/kZ6bjE++jMOa4/ewjg/N0zq0wpOlsbajkrUKCldLH3++ecAKu+Hb9iwocoiuhKJBK1atcKGDRvUn/ChnTt3YubMmQgICIBYLMbo0aMVBRwAlJWVIT4+HoWFhUpf08DAAOvXr8fs2bMhCAK8vLywevVqhIaGaqIJRERaZWUiwYyBXni1nwd+i0rDxtNJSMjMxzenkvDdmdsY3qVyMHh7J/aQE/2byrfhBg4ciL1798La2lpTmXQOb8MRkS6SywWcuJmJb04mIeL2P0MP+re1w7T+nujd2paDwalJU/b7WycmpWzsWCwRka67cjcHG08n4eC1dMgffit0cLbA1P6eCPJ14mBwapLUWizNmTMHy5Ytg6mpKebMmVPrsatXr1Y9rY5jsURETcXd7EJ8d6ZyMHhRWeVgcBcrY7zSpxWCe7rBzLBOq2QRNUpqLZYGDhyIffv2wcrKCgMHDqz5YiIRjh07VrfEOozFEhE1NQ8KSrEz4g62nvsbWfmlAABzI32E+LnjlT6t4GBhpOWERPXH23ANiMUSETVVxWUV+DUqFRtPJyHpXgEAwEBPhBFdXDC1vyfaOphrOSFR3WmsWNqxYwdGjRoFExOTeodsKlgsEVFTJ5cLOHYjExtPJSHy738Ggz/dzg5T+3vC35ODwUn3aKxYsrOzQ1FREYYPH46XX34ZgwcPrjKNQHPEYomImpOo5Af49nQSDsVIFYPBO7pYYGr/1gjq6Ah9DgYnHaGxYqm8vByHDh3CDz/8gN9++w0mJiYYM2YMQkJC0Lt373oH10UsloioObpzvwDfnbmNny7eRXGZHEDlYPApfT3w0lOuMOVgcGrkGmTMUmFhIfbt24ddu3bhyJEjaNmyJRITE+t6OZ3FYomImrMHBaX4/vwdbDv3N+4XVA4GtzDSx8u93DGpdyvYczA4NVINNsA7KysLu3fvxoYNGxAXF6fSQrpNBYslIqLKweB7L6di0+kkJGVVDgaX6IkxsmvlzOBtOBicGhmNFkuPepR27tyJo0ePwtXVFWPHjkVISIjG1odrzFgsERH9Qy4XcCQuAxtPJeHinQeK7c9422Nqf0/4edhwMDg1ChorloKDg3HgwAGYmJjgxRdfREhICPz9/esdWJexWCIierxLdx7g21NJOHxdikffNp1bWiK0vyeGdOBgcNIuZb+/VR59p6enh59++olPwRER0RN1d7dG9/HdcTurAN+dScLPF1NwJSUXM3dFwdXGGFP6eODFp1xhIuFgcGq8OCmlGrBniYhIOffzS/D9+TvYHn4H2Q8Hg1saG2B8L3dM7N0KduaGWk5IzYlab8OtXbsWU6dOhZGREdauXVvrsbNmzVI9rY5jsUREpJqi0gr8cjkFm04n4e/7hQAAib4Ygzs44lkfBwxoawdLYwMtp6SmTq3FkoeHBy5evAhbW1t4eHjUfDGRCElJSXVLrMNYLBER1U2FXEDY9QxsPJWIy8k5iu36YhF6etggsL0DAts7wM2Wq0aQ+nFtuAbEYomIqP6i7+bgcKwUR65n4FZmfpV97RzMEdDeHoE+DujS0gpiMZ+mo/rTWLG0dOlSvP3229XWhisqKsKqVauwePHiuiXWYSyWiIjU6++sAhyJy8CRuAxc+PsBKuT/fFW1MDNEgLc9Atrbo2+bFhwcTnWmsWJJT08P6enpsLe3r7L9/v37sLe356SULJaIiNQqt7AMJ25mIux6Bk7G30NeSblin6G+GH29WiDQxwEB3vacLZxUorGpAwRBeOxkYleuXIGNjY2qlyMiIqqVpYkBRnRxwYguLigtl+PC39kIu17Z65TyoAhHb2Ti6I1MAJVzOAW2d0CgjwO8Hc05+SWphdI9S9bW1hCJRIrq69+/gBUVFcjPz8drr72G9evXayxsY8WeJSKihicIAuIz8nA0rrLXKfpuTpX9LlbGCHw4zsnPwxYSfU6ASVWp/Tbctm3bIAgCJk+ejDVr1sDS0lKxTyKRoFWrVs12Jm8WS0RE2peZV4xjcZk4EpeJMwn3UFwmV+wzM9THgLZ2CPSxx8B29rAykWgxKTUWGhuzdPLkSfTu3RsGBpz/4hEWS0REjUtRaQXOJmThSFwGjt7IxL28EsU+PbEIPdytFbfrPFqYajEpaZNaiyWZTKa4iEwmq/XY5lgssFgiImq85HIBV1NzceThOKcb0rwq+1vbmSLQp3I+p25u1tDjtATNhlqLpX8/AScWix87YO7RwG8+DcdiiYioMbubXYijcRk4EpeJ80n3Uf6vaQlsTCUY2M4ez/rYo18bO5gaclqCpkytxdLJkyfRp08f6Ovr4+TJk7UeO2DAANXT6jgWS0REuklWXIZTN+/hyPUMHLuRCVnxP9MSSPTE8G9t+7DXyR5OlsZaTEqawBm8GxCLJSIi3VdWIcfFvx/gaFwGwuIycOfhmnWPdHC2QGB7Bzzr44AOzhaclqAJ0FixdOjQIZiZmaFv374AgPXr1+Pbb7+Fj48P1q9fD2tr6/ol10EsloiImhZBEJB4Lx9h1zNxJC4Dl5Mf4N/flo4WRgj0sUdAewf4e9rCyEBPe2GpzjRWLPn6+mLlypUICgrCtWvX0KNHD8ydOxfHjx+Ht7c3tmzZUu/wuobFEhFR05aVX4LjNyoLp1M3s1BU9s/4XBOJHvq3sUOgjwMGtrODrZmhFpOSKjRWLJmZmSEmJgatWrXChx9+iJiYGOzZsweXL19GUFAQpFJpvcPrGhZLRETNR3FZBcKT7iuersuQ/TMtgUgEdHezVjxd19rOlLfrGjGNLXcikUhQWFh5H/fIkSOYMGECAMDGxuaJ0woQERHpOiMDPQxsVzm55f9GdkRMqkyx6G9smgwX7zzAxTsP8PHBG2hla6KYz6mHuzX09TiLuC5SuWdp+PDhKC0tRZ8+fbBs2TLcvn0bLi4u+OuvvzBz5kzcvHlTU1kbLfYsERERAKTlFCmmJQhPvI/Sin9mETeV6KGDiyU6uVjCt6UlOrW0gruNCcSc10lrNHYbLjk5Ga+//jru3r2LWbNmYcqUKQCA2bNno6KiAmvXrq1fch3EYomIiP4rv6Qcp2/eQ1hcBo7fyMSDwrJqx5gb6cPXxbLy1dISnVys4GpjzFt3DYRTBzQgFktERFSbCnnl03VXU3JxLSUHV1NzcT1NhpJyebVjLY0N0KllZQHVqaUlfFtawdnSiAWUBmi0WJLL5UhISEBmZibk8qofdP/+/VVPq+NYLBERkarKKuS4lZGPa6k5lUVUai5upOdVuXX3iK2pBB0fFU8ulbfwHCwMWUDVk8aKpfPnz2PcuHG4c+cO/nsqlzthsURERHVXWi7HzYy8h8VTZREVL82rsiTLI3bmhv8a/2SJji6WsDc30kJq3aWxYqlLly5o27YtlixZAicnp2pVraWlZd0SP0F2djbeeOMN7N+/H2KxGKNHj8YXX3wBMzOzGs95+umnqy3PMm3aNGzYsEHxc3JyMqZPn47jx4/DzMwMEydOxIoVK6Cvr/yDgiyWiIhIU4rLKnBDmld5++5hD9TNjDw8pn6Co4XRw7FPlUWUr4sl532qhcamDrh16xb27NkDLy+vegVUVUhICNLT0xEWFoaysjK88sormDp1Knbt2lXreaGhoVi6dKniZxMTE8WfKyoqMGzYMDg6OuLcuXNIT0/HhAkTYGBggOXLl2usLURERMoyMtBDF1crdHG1UmwrKq3A9fRcXEvJxdXUyv9NuJcPqawY0uvFCLueoTjWxcr44dgnS8VgcisTiRZaortU7ll65pln8M4772DIkCGaylRNXFwcfHx8cOHCBfTo0QNA5bIrQUFBSElJgbOz82PPe/rpp9GlSxesWbPmsfsPHjyI5557DmlpaXBwcAAAbNiwAe+++y7u3bsHiUS5Xyb2LBERkbYVlJQjNk2Gqyk5uPawgErKKnjssW42JlV6oDq6WMLCyKCBE2ufxnqW3njjDcydOxdSqRS+vr4wMKj6l9upUyfV0z5BeHg4rKysFIUSAAQGBkIsFiMiIgIvvPBCjefu3LkTO3bsgKOjI55//nksWrRI0bsUHh4OX19fRaEEAIMHD8b06dMRGxuLrl27PvaaJSUlKCn5Z8ZWTsZJRETaZmqoj54eNujpYaPYJisuQ2yqrMog8jv3C5GcXfn642q64ljPFqZVBpF3cLGEmaHKZUKTpPLfwujRowEAkydPVmwTiUQQBEFjA7ylUins7e2rbNPX14eNjU2ty6uMGzcO7u7ucHZ2xtWrV/Huu+8iPj4ee/fuVVz334USAMXPtV13xYoVWLJkSV2bQ0RE1CAsjAzg39oW/q1tFdtyCksRkyrD1dScytt4KblIzSlCUlYBkrIK8PuVNACVS7e0tjOrMojcx8kSxpLmt2iwysXS7du31fbm8+fPx8qVK2s9Ji4urs7Xnzp1quLPvr6+cHJyQkBAABITE9G6des6X3fBggWYM2eO4meZTAZXV9c6X4+IiKihWJlI0LdNC/Rt00Kx7X5+Ca6l5iImNVfRA5WeW4yEzHwkZOZjb1QqAEAsAtrYmyuKJ18XS7R3soCRQdMuoFQultzd3dX25nPnzsWkSZNqPcbT0xOOjo7IzMyssr28vBzZ2dlwdHRU+v38/PwAAAkJCWjdujUcHR0RGRlZ5ZiMjMpBcbVd19DQEIaGfLqAiIiaBlszQzzdzh5Pt/vnLk5mXvE/xVNKLq6k5CIrvwTxGXmIz8jDnkspAAB9sQhtHMzhZW+G1namaG1nhtZ2ZvC0M20yRVSdbkZ+//332LBhA27fvo3w8HC4u7tjzZo18PDwwIgRI5S+jp2dHezs7J54nL+/P3JycnDp0iV0794dAHDs2DHI5XJFAaSM6OhoAICTk5Piuh999BEyMzMVt/nCwsJgYWEBHx8fpa9LRETU1NibG+EZbyM84105PEUQBGTISnA1JaeyiHpYSGUXlCIuXYa49Krjd0WiyifxHhVOj4qo1vamsDPTrQk1VX4a7uuvv8bixYvx1ltv4aOPPkJMTAw8PT2xdetWbNu2DcePH9dI0KFDhyIjIwMbNmxQTB3Qo0cPxdQBqampCAgIwPbt29GzZ08kJiZi165dCAoKgq2tLa5evYrZs2ejZcuWirmXKioq0KVLFzg7O+OTTz6BVCrF+PHj8eqrr6o0dQCfhiMiouZIEASk5RYjNrXyybvEzHwk3stH4r0C5BZVXwvvEXMj/WpFlJe9KdxsTCHRFzdYfo1NSunj44Ply5dj5MiRMDc3x5UrV+Dp6YmYmBg8/fTTyMrKqnf4x8nOzsbMmTOrTEq5du1axaSUf//9Nzw8PHD8+HE8/fTTuHv3Ll5++WXExMSgoKAArq6ueOGFF7Bw4cIqfyF37tzB9OnTceLECZiammLixIn4+OOPOSklERFRHQmCgPsFpUi6V1BZPD0sopKyCnA3u/CxE2oCgJ5YBHcbk2o9Ua3tzDQyN5TGiiVjY2PcuHED7u7uVYqlW7duoVOnTigqKqp3eF3DYomIiEg5xWUVuHO/sFoRlZiZj4LSmp+o3/t6b3Rzs1ZrFo3Ns+Th4YHo6OhqA70PHTqE9u3bq56UiIiImg0jAz20czRHO0fzKtsfjYmqvI2XX6VXKi23GB62plpKXIdiac6cOZgxYwaKi4shCAIiIyPxww8/YMWKFdi0aZMmMhIREVETJxKJ4GhpBEdLI/TxalFlX2FpOUwk2psgU+V3fvXVV2FsbIyFCxeisLAQ48aNg7OzM7744gsEBwdrIiMRERE1Y9oslIA6jFn6t8LCQuTn51ebXbu54ZglIiIi3aPs97fKz+cVFRWhsLAQAGBiYoKioiKsWbMGf/31V93TEhERETVSKhdLI0aMwPbt2wEAOTk56NmzJz777DOMGDECX3/9tdoDEhEREWmTysXS5cuX0a9fPwDAnj174OjoiDt37mD79u1Yu3at2gMSERERaZPKxVJhYSHMzSsf9/vrr78watQoiMVi9OrVC3fu3FF7QCIiIiJtUrlY8vLywq+//oq7d+/i8OHDGDRoEAAgMzOTg5uJiIioyVG5WFq8eDHefvtttGrVCn5+fvD39wdQ2cvUtWtXtQckIiIi0qY6TR0glUqRnp6Ozp07QyyurLciIyNhYWEBb29vtYds7Dh1ABERke7R2HInAODo6AhHR8cq23r27FmXSxERERE1airfhiMiIiJqTlgsEREREdWCxRIRERFRLVgsEREREdVCu8v4NhGPHiiUyWRaTkJERETKevS9/aSJAVgsqUFeXh4AwNXVVctJiIiISFV5eXmwtLSscX+d5lmiquRyOdLS0mBubg6RSKS268pkMri6uuLu3bvNbv6m5tr25tpugG1vjm1vru0G2PbG0nZBEJCXlwdnZ2fFvJGPw54lNRCLxWjZsqXGrm9hYaH1Xyhtaa5tb67tBtj25tj25tpugG1vDG2vrUfpEQ7wJiIiIqoFiyUiIiKiWrBYasQMDQ3xwQcfwNDQUNtRGlxzbXtzbTfAtjfHtjfXdgNsu661nQO8iYiIiGrBniUiIiKiWrBYIiIiIqoFiyUiIiKiWrBYIiIiIqoFi6VGbP369WjVqhWMjIzg5+eHyMhIbUeqlw8//BAikajKy9vbW7G/uLgYM2bMgK2tLczMzDB69GhkZGRUuUZycjKGDRsGExMT2NvbY968eSgvL2/optTq1KlTeP755+Hs7AyRSIRff/21yn5BELB48WI4OTnB2NgYgYGBuHXrVpVjsrOzERISAgsLC1hZWWHKlCnIz8+vcszVq1fRr18/GBkZwdXVFZ988ommm/ZET2r7pEmTqv0ODBkypMoxutj2FStW4KmnnoK5uTns7e0xcuRIxMfHVzlGXb/fJ06cQLdu3WBoaAgvLy9s3bpV082rlTJtf/rpp6t97q+99lqVY3Sx7V9//TU6deqkmFzR398fBw8eVOxvqp/5k9rdJD9vgRql3bt3CxKJRNi8ebMQGxsrhIaGClZWVkJGRoa2o9XZBx98IHTo0EFIT09XvO7du6fY/9prrwmurq7C0aNHhYsXLwq9evUSevfurdhfXl4udOzYUQgMDBSioqKEP//8U2jRooWwYMECbTSnRn/++afw/vvvC3v37hUACPv27auy/+OPPxYsLS2FX3/9Vbhy5YowfPhwwcPDQygqKlIcM2TIEKFz587C+fPnhdOnTwteXl7C2LFjFftzc3MFBwcHISQkRIiJiRF++OEHwdjYWPjmm28aqpmP9aS2T5w4URgyZEiV34Hs7Owqx+hi2wcPHixs2bJFiImJEaKjo4WgoCDBzc1NyM/PVxyjjt/vpKQkwcTERJgzZ45w/fp1Yd26dYKenp5w6NChBm3vvynT9gEDBgihoaFVPvfc3FzFfl1t+++//y788ccfws2bN4X4+HjhvffeEwwMDISYmBhBEJruZ/6kdjfFz5vFUiPVs2dPYcaMGYqfKyoqBGdnZ2HFihVaTFU/H3zwgdC5c+fH7svJyREMDAyEn3/+WbEtLi5OACCEh4cLglD5RSwWiwWpVKo45uuvvxYsLCyEkpISjWavq/8WDHK5XHB0dBRWrVql2JaTkyMYGhoKP/zwgyAIgnD9+nUBgHDhwgXFMQcPHhREIpGQmpoqCIIgfPXVV4K1tXWVdr/77rtCu3btNNwi5dVULI0YMaLGc5pK2zMzMwUAwsmTJwVBUN/v9zvvvCN06NChynu99NJLwuDBgzXdJKX9t+2CUPnl+eabb9Z4TlNpuyAIgrW1tbBp06Zm9ZkLwj/tFoSm+XnzNlwjVFpaikuXLiEwMFCxTSwWIzAwEOHh4VpMVn+3bt2Cs7MzPD09ERISguTkZADApUuXUFZWVqXN3t7ecHNzU7Q5PDwcvr6+cHBwUBwzePBgyGQyxMbGNmxD6uj27duQSqVV2mlpaQk/P78q7bSyskKPHj0UxwQGBkIsFiMiIkJxTP/+/SGRSBTHDB48GPHx8Xjw4EEDtaZuTpw4AXt7e7Rr1w7Tp0/H/fv3FfuaSttzc3MBADY2NgDU9/sdHh5e5RqPjmlM/y78t+2P7Ny5Ey1atEDHjh2xYMECFBYWKvY1hbZXVFRg9+7dKCgogL+/f7P5zP/b7kea2ufNhXQboaysLFRUVFT5RQIABwcH3LhxQ0up6s/Pzw9bt25Fu3btkJ6ejiVLlqBfv36IiYmBVCqFRCKBlZVVlXMcHBwglUoBAFKp9LF/J4/26YJHOR/Xjn+3097evsp+fX192NjYVDnGw8Oj2jUe7bO2ttZI/voaMmQIRo0aBQ8PDyQmJuK9997D0KFDER4eDj09vSbRdrlcjrfeegt9+vRBx44dFbnU8ftd0zEymQxFRUUwNjbWRJOU9ri2A8C4cePg7u4OZ2dnXL16Fe+++y7i4+Oxd+9eALrd9mvXrsHf3x/FxcUwMzPDvn374OPjg+jo6Cb9mdfUbqBpft4slqjBDB06VPHnTp06wc/PD+7u7vjpp5+0/o88NYzg4GDFn319fdGpUye0bt0aJ06cQEBAgBaTqc+MGTMQExODM2fOaDtKg6up7VOnTlX82dfXF05OTggICEBiYiJat27d0DHVql27doiOjkZubi727NmDiRMn4uTJk9qOpXE1tdvHx6dJft68DdcItWjRAnp6etWemsjIyICjo6OWUqmflZUV2rZti4SEBDg6OqK0tBQ5OTlVjvl3mx0dHR/7d/Jony54lLO2z9bR0RGZmZlV9peXlyM7O7tJ/V0AgKenJ1q0aIGEhAQAut/2mTNn4sCBAzh+/Dhatmyp2K6u3++ajrGwsND6f3DU1PbH8fPzA4Aqn7uutl0ikcDLywvdu3fHihUr0LlzZ3zxxRdN/jOvqd2P0xQ+bxZLjZBEIkH37t1x9OhRxTa5XI6jR49WuSes6/Lz85GYmAgnJyd0794dBgYGVdocHx+P5ORkRZv9/f1x7dq1Kl+mYWFhsLCwUHT/NnYeHh5wdHSs0k6ZTIaIiIgq7czJycGlS5cUxxw7dgxyuVzxj46/vz9OnTqFsrIyxTFhYWFo166d1m9DqSIlJQX379+Hk5MTAN1tuyAImDlzJvbt24djx45Vu02ort9vf3//Ktd4dIw2/114UtsfJzo6GgCqfO662PbHkcvlKCkpadKf+eM8avfjNInPWyvDyumJdu/eLRgaGgpbt24Vrl+/LkydOlWwsrKq8vSArpk7d65w4sQJ4fbt28LZs2eFwMBAoUWLFkJmZqYgCJWP2bq5uQnHjh0TLl68KPj7+wv+/v6K8x89bjpo0CAhOjpaOHTokGBnZ9fopg7Iy8sToqKihKioKAGAsHr1aiEqKkq4c+eOIAiVUwdYWVkJv/32m3D16lVhxIgRj506oGvXrkJERIRw5swZoU2bNlUen8/JyREcHByE8ePHCzExMcLu3bsFExMTrU8dUFvb8/LyhLffflsIDw8Xbt++LRw5ckTo1q2b0KZNG6G4uFhxDV1s+/Tp0wVLS0vhxIkTVR6XLiwsVByjjt/vR49Tz5s3T4iLixPWr1+v9cfIn9T2hIQEYenSpcLFixeF27dvC7/99pvg6ekp9O/fX3ENXW37/PnzhZMnTwq3b98Wrl69KsyfP18QiUTCX3/9JQhC0/3Ma2t3U/28WSw1YuvWrRPc3NwEiUQi9OzZUzh//ry2I9XLSy+9JDg5OQkSiURwcXERXnrpJSEhIUGxv6ioSHj99dcFa2trwcTERHjhhReE9PT0Ktf4+++/haFDhwrGxsZCixYthLlz5wplZWUN3ZRaHT9+XABQ7TVx4kRBECqnD1i0aJHg4OAgGBoaCgEBAUJ8fHyVa9y/f18YO3asYGZmJlhYWAivvPKKkJeXV+WYK1euCH379hUMDQ0FFxcX4eOPP26oJtaotrYXFhYKgwYNEuzs7AQDAwPB3d1dCA0NrfYfALrY9se1GYCwZcsWxTHq+v0+fvy40KVLF0EikQienp5V3kMbntT25ORkoX///oKNjY1gaGgoeHl5CfPmzasy744g6GbbJ0+eLLi7uwsSiUSws7MTAgICFIWSIDTdz7y2djfVz1skCILQcP1YRERERLqFY5aIiIiIasFiiYiIiKgWLJaIiIiIasFiiYiIiKgWLJaIiIiIasFiiYiIiKgWLJaIiIiIasFiiYiIiKgWLJaIiIiIasFiiYioFrNnz8aoUaO0HYOItIjFEhFRLSIjI9GjRw9txyAiLeLacEREj1FaWgpTU1OUl5crtvn5+eH8+fNaTEVE2qCv7QBERI2Rvr4+zp49Cz8/P0RHR8PBwQFGRkbajkVEWsBiiYjoMcRiMdLS0mBra4vOnTtrOw4RaRHHLBER1SAqKoqFEhGxWCIiqkl0dDSLJSJisUREVJNr166hS5cu2o5BRFrGYomIqAZyuRzx8fFIS0tDbm6utuMQkZawWCIiqsH//vc/bN26FS4uLvjf//6n7ThEpCWcZ4mIiIioFuxZIiIiIqoFiyUiIiKiWrBYIiIiIqoFiyUiIiKiWrBYIiIiIqoFiyUiIiKiWrBYIiIiIqoFiyUiIiKiWrBYIiIiIqoFiyUiIiKiWrBYIiIiIqrF/wNHTA6OTX7COQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "solution = sim.solve(\n",
    "    [0, 3600], inputs={\"Current function [A]\": 0.15652}, calculate_sensitivities=True\n",
    ")\n",
    "plt.plot(\n",
    "    solution.t, solution[\"Terminal voltage [V]\"].sensitivities[\"Current function [A]\"]\n",
    ")\n",
    "\n",
    "plt.xlabel(r\"$t$\")\n",
    "plt.ylabel(\"sensitivities of Terminal voltage wrt Current fuction\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d84c83fb",
   "metadata": {},
   "source": [
    "## Sensitivities and data fitting\n",
    "\n",
    "Sensitivities are often used to aid data fitting by providing a means to calculate the gradient of the function to be minimised. Take for example the data fitting exercise we introduced in the previous notebook. Once again we will generate some fake data for fitting, like so:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "33d95c39",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_eval = np.linspace(0, 3600, 100)\n",
    "data = sim.solve([0, 3600], inputs={\"Current function [A]\": 0.2222})[\n",
    "    \"Terminal voltage [V]\"\n",
    "](t_eval)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5045b57",
   "metadata": {},
   "source": [
    "Now we will contruct a function to minimise, but here we will return both the value of the function, and its gradient with respect to the input parameter \"Current function\". Note that our objective function is the sum of squared different between the vector $\\mathbf{f}$, the simulated \"Terminal voltage\", and $\\mathbf{d}$, the vector of fake data, given by\n",
    "\n",
    "$$\\mathcal{O}(a) = \\sum_{i=0}^{i=N} (f_i(a) - d_i)^2$$ \n",
    "\n",
    "where $a$ is the parameter to be optimised (in this case \"Current function\"), $f_i$ is each element of the vector $\\mathbf{f}$, and $d_i$ is each element of $\\mathbf{d}$. We wish to also find the gradient of this function wrt the parameter $a$, which is:\n",
    "\n",
    "$$\\frac{\\partial \\mathcal{O}}{\\partial a}(a) = 2 \\sum_{i=0}^{i=N} (f_i(a) - d_i) \\frac{\\partial f_i}{\\partial a} $$ \n",
    "\n",
    "Using these equations, we will define a function that takes in as an argument $a$, and returns $(\\mathcal{O}(a), \\frac{\\partial \\mathcal{O}}{\\partial a}(a))$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ffad2bc0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sum_of_squares_jac(parameters):\n",
    "    sol = sim.solve(\n",
    "        [0, 3600],\n",
    "        t_interp=t_eval,\n",
    "        inputs={\"Current function [A]\": parameters[0]},\n",
    "        calculate_sensitivities=True,\n",
    "    )\n",
    "    term_voltage = sol[\"Terminal voltage [V]\"].data\n",
    "    term_voltage_sens = sol[\"Terminal voltage [V]\"].sensitivities[\n",
    "        \"Current function [A]\"\n",
    "    ]\n",
    "\n",
    "    f = np.sum((term_voltage - data) ** 2)\n",
    "    g = 2 * np.sum((term_voltage - data) * term_voltage_sens)\n",
    "    print(\n",
    "        f\"evaluating function and jacobian for p = {parameters[0]}, \\tloss = {f}, grad = {g}\"\n",
    "    )\n",
    "    return f, g"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fdcae8ac",
   "metadata": {},
   "source": [
    "We can then use this function along with an optimisation algorithm to recover the value of the Current function that was used to generate the data. In this case we will use the `scipy.optimize` module again. This module allows the use of a function in the form given above to perform the minimisation, using both the value of the objective function and its gradient to find the minimum value of $a$ in the least number of steps.\n",
    "\n",
    "Once again, we will place bounds on \"Current function [A]\" between $(0.01, 0.6)$, and use a random starting value $x_0$ between these bounds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "44f52a7e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "starting parameter is 0.4035076613514513\n",
      "evaluating function and jacobian for p = 0.4035076613514513, \tloss = 0.358632833514908, grad = 3.7278265436340283\n",
      "evaluating function and jacobian for p = 0.01, \tloss = 0.8765036816401419, grad = -9.691403058800336\n",
      "evaluating function and jacobian for p = 0.23970725060294026, \tloss = 0.0036706826105448887, grad = 0.41303112361463107\n",
      "evaluating function and jacobian for p = 0.21929734316448973, \tloss = 0.00010628748342624681, grad = -0.07293742235816741\n",
      "evaluating function and jacobian for p = 0.22236059911277223, \tloss = 2.5627985467376454e-07, grad = 0.0035066000930420284\n",
      "evaluating function and jacobian for p = 0.22222008304298907, \tloss = 7.758859239380127e-09, grad = 2.935984691530423e-05\n",
      "evaluating function and jacobian for p = 0.22221889660489277, \tloss = 7.74141508243804e-09, grad = -1.183255640750795e-08\n",
      "recovered parameter is 0.22221889660489277\n"
     ]
    }
   ],
   "source": [
    "bounds = (0.01, 0.6)\n",
    "x0 = np.random.uniform(low=bounds[0], high=bounds[1])\n",
    "\n",
    "print(\"starting parameter is\", x0)\n",
    "res = scipy.optimize.minimize(sum_of_squares_jac, [x0], bounds=[bounds], jac=True)\n",
    "print(\"recovered parameter is\", res.x[0])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "env",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
